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Abstract 



A previously-proposed method of constructing spatially-extended gauge-invariant 
three-quark operators for use in Monte Carlo lattice QCD calculations is tested, and a 
methodology for using these operators to extract the energies of a large number of baryon 
states is developed. This work is part of a long-term project undertaken by the Lattice 
Hadron Physics Collaboration to carry out a first-principles calculation of the low-lying 
spectrum of QCD. 

The operators are assemblages of smeared and gauge-covariantly-displaced quark fields 
having a definite flavor structure. The importance of using smeared fields is dramatically 
demonstrated. It is found that quark field smearing greatly reduces the couplings to the 
unwanted high-lying short-wavelength modes, while gauge field smearing drastically 
reduces the statistical noise in the extended operators. Group-theoretical projections onto 
the irreducible representations of the symmetry group of a cubic spatial lattice are used to 
endow the operators with lattice spin and parity quantum numbers, facilitating the 
identification of the J p quantum numbers of the corresponding continuum states. 

The number of resulting operators is very large; consequently a key aspect of this work is 
the development of a selection method for finding a sufficient subset of operators for 
accurately extracting the lowest seven or eight energy levels in each symmetry channel. A 
procedure in which the diagonal elements of the correlation matrix of the operators are 
first evaluated to remove noisy operators, followed by the selection of sixteen operators 
whose renormalized correlation matrix at a fixed small time separation has a low condition 
number for both the even- and odd-parity channels, is found to work well. 

These techniques are applied in the construction of nucleon operators. Correlation matrix 
elements between these operators are estimated using 200 configurations on a 12 3 x 48 
anisotropic lattice in the quenched approximation with unphysically heavy u, d quark 
masses (the pion mass is approximately 700 MeV). After a change of basis operators using 
a variational method is applied, the energies of up to eight states are extracted in each 
symmetry channel. Although comparison with experiment is not justified, the pattern of 
levels obtained qualitatively agrees with the observed spectrum. A comparison with quark 
model predictions is also made; the quark model predicts more low-lying even-parity states 
than this study yields, but both the quark model and this study predict more odd-parity 
states near 2 GeV than currently observed in experiments. 
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Chapter 1 
Introduction 



Particle physicists seek to identify the elementary building blocks of the universe, and the 
mechanisms by which these elements interact. This view, known as reductionism, is seen 
by many as the necessary starting point for the discovery of a unified theory which 
quantitatively describes all natural phenomena. If we can identify the building blocks and 
the rules which hold them together, then we can begin to explore which features of our 
universe arise from the complex interactions of these fundamental components. 



1.1 The atom and nucleus 



In 1803, John Dalton presented his atomic theory which stated that all compounds were 
composed of and reducible to collections of atoms. This sparked a great effort to use any 
chemical means available to separate compounds into their fundamental atoms, the so 
called elements. As scientists discovered new elements and observed their properties in 
reactions, they began to notice that entire groups of elements behaved in similar ways. In 
1869, Dmitri Mendeleev |HM84] introduced his now ubiquitous periodic table of the 
elements, which placed elements with similar reaction properties into columns, ordered by 
atomic mass. Mendeleev's table classified the 63 known elements at the time, and he 
predicted the existence and characteristics of three elements which were later discovered: 
gallium, scandium, and germanium. 

The similar behavior of the different elements in each column was a reassurance that there 
was indeed an order hidden beneath the diversity of the everyday world. Today, the 118 
known elements fit neatly into the eighteen columns of the modern periodic table (there are 
also the lanthanoid and actinoid classifications). These elements include stable atoms such 
as hydrogen, carbon, nitrogen, and oxygen; and unstable atoms, such as radium and 
uranium which last long enough to be characterized, but which eventually break down into 
lighter elements through radioactive decay. 
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The large size and redundancy in the table were also significant hints that the atom 
possessed substructure. A breakthrough on this front was the identification of the electron 
in 1897 by J.J. Thomson |Gri87] . At this time, scientists were fascinated by the cathode 
ray tube, a glass chamber with a metal electrode at each end and a port for a vacuum 
pump. When the air was pumped out of the tube and a voltage difference was applied to 
the electrodes, the tube would glow with 'rays' emanating from the negative electrode (the 
cathode). This effect was later enhanced by J.B. Johnson who added a heating element to 
the cathode. Before 1897 scientists suspected, but were not certain, that the rays were 
composed of electrically charged particles. 

J.J. Thomson's cathode ray experiments built on those of his contemporaries, and involved 
the deflection of these rays by known electric and magnetic fields. By balancing the electric 
force Fe = qE with the magnetic force Fb = qv x B, he determined the velocity of the 
particles to be roughly one tenth of the speed of light from the relation 

E 



He then switched off the electric field, measured the radius of curvature R of the beam, and 
used the condition of uniform circular motion F = mv 2 /R to determine the charge-to-mass 
ratio from the relation 

q v 
m ~ KB' 

Thomson found that the charge-to-mass ratio was abnormally large when compared to the 
known ions. He assumed that the mass was very small (as opposed to assuming that the 
charge was very large), and concluded that the ray was composed of negatively charged 
elementary particles, and identified the particles as constituents of the atom. Later, they 
were named 'electrons^, taken from a term coined in 1894 by electro-chemist G. Johnstone 
Stoney. Thomson had discovered that cathode ray electrons came from the atoms 
composing the cathode; these electrons were 'boiled' off of the cathode as it heated up. 

The charge of the electron (and mass, from Thomson's ratio) was measured in 1909 by 
Robert Millikan in his oil-drop experiment. He used an atomizer to spray a mist of oil 
droplets above two plates. The top plate had a small hole through which a few droplets 
could pass. Millikan would vary the electric field between the plates until he had one 
(charged) drop of oil suspended. Then he turned off the electric field and allowed the drop 
to fall until it had reached its terminal velocity V\. The velocity-dependent drag force on 
the drop was given by Stoke's Law: 

F D1 = 6-Krrjvx, 

where r is the radius of the (assumed spherical) drop and rj is the viscosity of air. At 
terminal velocity, the drag force balanced the the 'apparent weight' of the oil drop 

1 From the Greek elektron, meaning amber. Early experiments in electricity used amber rods rubbed with 
fur to build up a static charge. 
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W, which is 



W= (4/3)7rr 3 (p oil -p air )(7, 



where p ii is the density of the oil, and p a ; r is the density of air. His measurement of v\ 
allowed him to infer r, and thus the apparent weight of the drop W. He then turned the 
electric field back on before the drop reached the bottom plate and measured the terminal 
velocity v 2 of the drop's upward motion. At this point, Millikan knew that 
qE = F D2 + W = F m v 2 /vi + W and could therefore determine the charge q on the drop: 



After repeating the experiment for many droplets, Millikan confirmed that the quantity of 
charge on a drop was always a multiple of the same number, the charge on a single electron 
(1.602 x KT 19 Coulomb, in SI units). 

Electrically neutral atoms must possess enough positive charge to compensate for the 
negative charge of the electrons. The modern model of the atom was born from Geiger and 
Marsden's |GM09] alpha scattering experiments performed under Ernest 
Rutherford [Rut EL], a former student of Thomson's, in 1909. The goal of these experiments 
was to probe the positive charge distribution of the atom. Charged alpha particles (doubly 
ionized helium atoms) emitted from a radioactive radium source were directed at a gold 
foil. A zinc sulfide screen was placed at various positions to detect the scattered alpha 
particles. Rutherford, Geiger, and Marsden found that most of the alpha particles passed 
through the foil with little deflection but some deflected through large angles. This 
suggested that the positive charge in each atom was concentrated at the center and 
occupied just a fraction of the total atomic volume. 

The nucleus of the lightest element (Hydrogen) was named the protorJl. Its charge is 
exactly equal to the magnitude of the electron charge, but its mass is roughly 2000 times 
greater. In the atomic model, a hydrogen atom consisted of a proton and electron bound 
together via the electromagnetic force. 



1.2 Hydrogen spectroscopy 

Once scientists had a model of the hydrogen atom, they were in a position to discuss the 
hydrogen spectrum. When an electrical current is passed through pure hydrogen gas, the 
atoms absorb energy and then radiate at specific discrete wavelengths, which can be 
observed by passing the emitted light through a prism or diffraction grating. Spectroscopy 
was pioneered by people such as A.J. Angstrom and was in widespread use in chemistry for 
the classification of elements. The first observed series of emission lines was the Balmer 
series, named after J.J. Balmer who in 1885 first developed the empirical relationship for 




2 From the Greek proton, meaning first. 
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the spacing of the lines. The relationship was generalized in 1888 by Johannes Rydberg for 
the complete emission line spectrum: 



where k = 1/A is the discrete emission wavenumber and 1Z = 10967757.6 ± 1.2 m 1 is the 
Rydberg constant for the hydrogen spectrum [ ER85j . The Balmer series lines correspond to 



A burning question was: could the atomic model reproduce the observed hydrogen 
spectrum? 

In 1913, Niels Bohr combined the Rutherford model of the atom with Einstein's quantum 
theory of the photoelectric effectS introduced in 1905 by introducing his quantization 
condition. Bohr's model related the electromagnetic radiation emitted by an atom to 
electron transitions from states of definite angular momentum: 



The quantization of atomic energy states was experimentally verified in 1914 by Frank and 
Hertz who accelerated electrons through a potential difference in a tube filled with mercury 
vapor [ER85J. They measured the current as a function of applied voltage, which was an 
indicator of the number of electrons which passed unimpeded through the gas. When the 
kinetic energy of the electrons reached a threshold level, the current abruptly dropped. 
Frank and Hertz interpreted this to mean that the electrons were exciting the mercury 
atoms and being scattered. As the voltage was increased, the current would increase until 
successive thresholds were reached. This showed that the mercury atoms possessed discrete 
energy levels. The mercury atoms would only absorb energy from the electrons by 
transitioning from one energy level to another. 

The urge to understand the underpinnings of Bohr's (and in 1916 Sommerfeld and 
Wilson's) quantization conditions sparked the rapid development of an entirely new field of 
physics: quantum mechanics. Two notable contributions were Werner Heisenberg's 1925 
matrix mechanics paper, and Erwin Schrodinger's 1926 paper on quantization as an 
eigenvalue problem (which introduced 'Schrodinger's equation'). Quantum mechanics 
succeeded not only in accurately describing atomic spectra, but also lead to revolutionary 
technological advances such as the semiconductor used in computers and nuclear magnetic 
resonance (NMR) used in medical imaging. 

The hydrogen spectrum not only provided a test of the atomic model, it also led to 
refinements in our understanding of the atom and opened up an entirely new branch of 
physics. The importance of spectroscopy cannot be overstated. 

3 Einstein proposed that radiant energy comes in quanta known as photons with the energy frequency 
relation E = hv. 




n% = 2 and 712 — 3, 4, 



L = nh. 
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1.3 Nucleons and the strong nuclear force 

The neutron, the electrically neutral partner of the proton in the nucleus, was identified in 
1932 by James Chadwick |Cha32j from previous experiments in which a polonium source 
was used to bombard beryllium with alpha particles. 

In modern notation |CG89] , the nuclear reaction was: 

H<4 + Be2-Cj a + nS, 

where the subscript denotes the atomic number (number of protons), and the superscript 
denotes the atomic weight (number of protons+neutrons) . 

The observation that every atom contains a nucleus in which protons and neutrons are 
confined within 0.01% of the volume of the atom raises the question: what keeps the 
protons from flying apart from electrostatic repulsion? Physicists inferred the existence of a 
new force which could overpower the electromagnetic force, but which had a range on the 
order of the nuclear radius. In 1934, Hideki Yukawa worked out the quantization of the 
strong nuclear force field and predicted a new particle, the pion [Gri87j . In Yukawa's 
theory, nucleons interact with each other via the pion field. When the field is quantized 
according to the formalism of quantum field theory, the potential felt by the nucleons goes 
as 

1 r m^r "i 
r \ he J 

where r is the inter-nucleon separation, and m n is the mass of the pion, the quantum which 
mediates the strong nuclear forc^l To get a force with a range of 1 fm, the order of the 
typical nuclear radius 

7?V « ~ 200 MeV. 
1 tm 

It was then up to particle physicists to use all of the tools at their disposal to confirm the 
existence of the pion and complete the atomic model. 



1.4 Particle sources, particle detectors, and the 
particle zoo 

1.4.1 Sources of subatomic particles 

We have already discussed the cathode ray tube, a source of electrons when the cathode is 
heated enough to boil them off of the metal. A wide variety of experiments were performed 
which involved accelerating the electrons through a potential difference. 

4 Note that the electromagnetic force is mediated by the massless photon, giving the usual 1/r Coulomb 
potential. 
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In 1895, Henri Becquerel discovered that uranium was radioactive from the darkening of 
photographic film [CG89] . The radioactive decay of heavy elements such as polonium 
produced light nuclei such as that of helium (two protons and two neutrons), the so called 
alpha radiation. When a neutron decays into a proton, it emits an electron, the so called 
beta radiation. Electromagnetic radiation is called gamma radiation. Sources of radiation 
would be placed in front of collimators, which would only let through thin beams of the 
particles. These beams could then be directed onto various targets and the reactions 
observed. 

In 1912, Victor Hess used a hot air balloon to take measurements of ionizing radiation at 
varying altitudes. He discovered that the rate of ionization was roughly four times greater 
at an altitude of 5,300 meters than it was at ground level, thus showing that the radiation 
which ionizes the atmosphere is cosmic in origin. The discovery of these so called cosmic 
rays led to the discovery of a wide range of subatomic particles. 

Particle accelerators ('atom smashers') use electromagnetic fields to accelerate beams of 
charged particles to velocities comparable with the speed of light. Accelerators may be 
linear or circulajf). Circular colliders use magnetic fields to curve the paths of particles as 
they are accelerated around the ring. These high energy beams of particles can be directed 
upon stationary targets or made to collide with other beams of particles. Particles of one 
type could also be directed onto targets, producing particles of a different type which can 
themselves be collimated into a beam and focused with magnets. 

1.4.2 Detectors of subatomic particles 

A conspicuous signature of a charged cosmic ray is the ion trail left behind as electrons are 
stripped off of atoms in the ray's path. Some particle detectors turn these ion trails into 
visible paths, but this is not the only way to detect a subatomic particle. Here is a brief 
description of some of the tools experimentalists use to detect subatomic particles. 

Nuclear emulsions are a special mixture of gelatin and silver bromide salts and work 
under the same principles as chemical photography. The track of a charged particle creates 
sites of silver atoms on the salt grains. Photographic development chemicals reduce the 
silver bromide salts to silver, but are most effective when there are already silver atoms 
present. Early cosmic ray researchers would stack plates of nuclear emulsions, and then 
develop them into pictures of subatomic paths. 

Cloud chambers, also called Wilson chambers after their inventor C.T.R. Wilson, contain 
a supersaturated vapor (such as isopropyl alcohol) that forms droplets along the trail of 
ionization. The trails are illuminated with a bright light source and photographed. 

Bubble chambers work under the same principle as cloud chambers, but use a 



5 The Large Hadron Collider (LHC) is a proton-proton collider with a ring 27 km in circumference. 
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superheated transparent liquid (such as liquid hydrogen). When a charged particle passes 
through the liquid, it interacts with the molecules and deposits enough energy to boil the 
liquid around the interaction points. The result is the formation of a string of small 
bubbles along the particle's path which can then be illuminated and photographed. 

Drift chambers are described by [Col02j: 

Charged particles passing through the chambers ionize the gas in the chambers 
(DME or an argon/ethane mixture), and the resulting ionization is collected by 
wires maintained at high voltage relative to the surrounding "field" wires. The 
electrical signals from these wires are amplified, digitized, and fed into a 
computer which reconstructs the path of the original particle from the wire 
positions and signal delay times. 

Silicon detectors operate under the same principle as drift chambers but utilize a 
semiconducting material instead of a gas. This allows for much a higher energy resolution 
and a much higher spatial resolution. However, they are more expensive and more sensitive 
to the degradating effects of radiation than drift chambers. 

Scintillators are compounds which absorb energy from interactions with charged 
particles, and then re-emit that energy in the form of electromagnetic radiation at a longer 
wavelength (fluorescence). Scintillators have short decay times and are optically 
transparent to the flashes. Thus, particle tracks register as a series of rapid flashes within 
the material which can be seen by light detectors (such as photomultiplier tubes). 

Calorimeters measure the energy content of the particle shower which occurs when a 
subatomic particle strikes a dense barrier in the detector. Often calorimeters are segmented 
into different chambers, and the energy deposited by the particle showers in each chamber 
is used to infer the particle's identity and direction of travel. Calorimeters can be designed 
to detect either electromagnetic showers or hadron showers, which occur when a subatomic 
particle interacts with the barrier. 

1.4.3 The particle zoo 

Particle physicists use detectors to examine the contents of cosmic rays and the products of 
beam-target and beam-beam scattering experiments. They measure the trajectories either 
directly in the detector, or through reconstruction via the principle of conservation of mass 
and momentum. A particularly useful quantity measured is the differential cross section, 
the reaction rate per unit incident flux as a function of angle, energy, and any other 
parameters of interest [Per 00] . Other properties of the particles such as mass, spin, parity, 
form factors (describing the structure of the particle), life-time, and branching ratios 
(describing the relative likelihood of decay into each of several final states), can then be 
deduced from interaction cross section data. 
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There are three classes of outcomes when two subatomic particles A and B approach each 
other: 

1. Nothing happens: the particles pass right by each other without interacting 

2. Elastic scattering (A + B —>■ A + B): the particles interact through the exchange of a 
'messenger' particle. The outgoing particles are of the same type as the incoming 
particles, but may have a different energy and trajectory. 

3. Inelastic scattering (A + B — > C + D + • • •): Einstein's mass-energy relational 

E = / ymc 2 allows for the transmutation the combined energy and mass of the input 
particles into a completely different set of output particles, subject to quantum 
transition rules and kinematic constraints. 

It was expected that the above sources and detectors of particles would lead to the 
discovery of Yukawa's pion, and the completion of the atomic model. Some scientists 
anticipated that the proton, neutron, electron, photon, and pion would constitute the 
fundamental building blocks of all matter. It was a great surprise, then, when detailed 
observations of cosmic rays and scattering experiments uncovered a plethora of subatomic 
particles. One of the early cosmic ray candidates for the pion ended up being the muon, a 
more massive relative of the electron. I.I. Rabi put it best when he said "Who ordered 
that?" fCG89] 

The pion was experimentally confirmed in 1947 when D.H. Perkins |Per47] observed the 
explosion of a nucleus after capturing a cosmic-ray pion. He saw the ion trails created by 
the incoming pion and outgoing nuclear debris in a photographic emulsion. Nuclear 
disintegration by pion capture had been predicted in 1940 by Tomonaga and Araki. 

Modern estimates [H + 02] of the mass of the charged and neutral pions are: 



Pion 


Mass (MeV) 




139.57018(35) 


7T° 


134.9766(6) 



which is close to the very rough estimate of 200 MeV made by Yukawa. 

The muon, now identified as a distinct particle from the pion, was only the beginning of a 
long revolutionary series of discoveries. Positrons, kaons, antiprotons, tau leptons, 
neutrinos, and more took their place with the proton, neutron, and electron in the rapidly 
expanding 'particle zoo.' 

Most of these particles were unstable resonances. In analogy with radioactive elements in 
the periodic table which decay into simpler elements, a resonance is defined as an object of 

6 Einstein's original formulation of special relativity made a distinction between the mass of a body at 
rest and the perceived mass of a body in motion. Modern convention takes mass to be a characteristic 
of an object (e.g. on equal footing with charge), and explicitly includes the relativistic dilatation factor 
7 = l/ v /l-i>/c 2 . 
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mass M with a lifetime r much longer than the period associated with its 'characteristic 
frequency' v = M/h: 

r » h/M. 

For very massive resonances, this lifetime can be very short and is quoted in terms of the 
decay width T = h/r, a. natural spread in the energy of the decaying state induced by the 
uncertainty principle |Per00j . Resonances are then defined by the property (neglecting 
factors of 2ir and working in energy units) 

T « M. 

For example the Z° resonance, discovered in 1983 at CERN, has a mass of approximately 
91 GeV [H + 02j and a decay width of approximately 2.5 GeV, which corresponds to a 
lifetime of roughly 2.6 x 10~ 25 seconds, enough time for light to travel about one-tenth of a 
Fermi, much less than than the spatial extent of a proton. Such resonances cannot, 
therefore, be observed directly in particle detectors. Rather, resonances in the elastic 
channeljj 

A + B -> X -> A + B, 

show up as enhancements in the differential and total cross sections compared to what 
would be expected from simple kinematics alone. Resonances in the inelastic channel 

A + B ^ X ^ C + D + ... 

can be identified from trajectory reconstruction and decay products. 

The overwhelming number of subatomic particles discovered (including resonances) lead 
Wolfgang Pauli to exclaim in the 1950s "Had I foreseen this, I would have gone into 
botany." 



1.5 A new periodic table 
1.5.1 Classification of particles 

As new particles were discovered, they were classified by their values of conserved 
quantities, such as mass and electric charge. The existence of 'forbidden' reactions, which 
were never observed, led to the introduction of two more conserved quantum numbers: 
lepton number and baryon number. Leptons (such as the electron and muon) have lepton 
number = 1, baryon number = 0, and do not participate in any strong interactions. 
Hadrons have lepton number = 0, participate in strong interactions, and either have 
baryon number = 1 for baryons (such as the proton and neutron), or baryon number = 

A channel refers to the way a reaction can proceed. The channel concept is also applied to initial, 
intermediate, or final states in a reaction, e.g. the 'nucleon channel' contribution for a process has a nucleon 
resonance as an intermediate state, and the A ++ can be created in the p-zr channel. 
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for mesongj (such as the pion). Because an antiparticle annihilates its corresponding 
particle, it must have the opposite quantum numbers (except for mass, which is always 
positive and which is conserved in the radiation emitted from an annihilation reaction). 

The discovery of kaons, unstable mesons which were created easily via the strong 
interaction but which decayed slowly via a different process (the weak interaction, which 
also governs radioactive beta decay), lead to the introduction of 'strangeness', a new 
quantum number. The strong interaction conserves strangeness, but the weak interaction 
violates it. The time-scale over which the strong interaction acts is very short, thus giving 
rapid creation and decay rates. On the other hand, weak interactions have long time-scales 
which cause strangeness-changing reactions to proceed slowly. For example, a kaon (meson 
with strangeness = +1) and lambda (baryon with strangeness = -1) can be produced easily 
by the strong reaction: 

tT +p + ^K° + A°, \as\=q, 
but must decay separately via the much slower weak reactions 

K° -> 7T + 7T-, \AS\ = 1, 

A -> p + 7T- } \AS\ = I. 

In the tradition of Mendeleev, the discovered hadrons were placed into tables according to 
their masses and reaction properties (quantum numbers) independently by Murray 
Gell-Mann and Yuval Ne'eman in 1961 [GMNO Oj. In this arrangement scheme, dubbed the 
'Eightfold Way' by Gell-Mann, particles with similar masses were placed into hexagonal 
and triangular arrays which were labeled by strangeness S and electric charge Q. Two 
examples of these tables are shown in Figure 11.11 

Just as Mendeleev had predicted the existence and properties of gallium, scandium, and 
germanium from holes in his periodic table, so also did Murray Gell-Man predict the 
existence, mass, and quantum numbers of a new particle, the Q~ , a strangeness = -3 heavy 
baryon, which was later discovered at Brookhaven National Laboratory in 1964 |B + 64| . 

The existence of such a large number of hadrons along with their classification according to 
a comparably small number of quantum numbers hinted, once again, that there was 
substructure yet to be discovered. What was needed was a model of nucleon substructure 
which could explain the patterns in Gell-Man and Ne'eman's tables. 

1.5.2 The constituent quark model 

The constituent quark model describes hadrons as being built up from combinations of 
point-like particles. This model was proposed independently in 1964 by Gell-Mann and 
Zweig, and Gell-Mann named the constituents 'quarks^. 

8 From the Greek mesos, meaning middle. Mesons are so named because the first mesons had masses 
between that of the electron (0.5 MeV) and the proton (938 MeV). 

9 From James Joyce's Finnigan's Wake, referring to the sound a seagull makes. 
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Quarks are spin- 1/2 objects (fermions) with baryon number = 1/3 and come in different 
flavors. Originally three flavors were introduced: 'up', 'down', and 'strange': 





Charge 


Strangeness 


u 


+2/3 





d 


-1/3 





s 


-1/3 


-1 



The quarks combine into baryon and meson multiplets following the well-established rules 
of addition of angular momenta. According to the quark model, baryons are three-quark 
states, and mesons are quark- ant iquark states. All of the hadrons known in 1964 (and most 
discovered since) could be described using this model. Figure 11.21 shows baryons and 
mesons as viewed in the quark model, and figure 11.31 shows pictorially the reaction 

7i- +p + -> K° + A . 



Hadron 


Quark Content 


Baryon Number 


Charge 


Strangeness 


p+ 


uud 


1 


+ 1 





n° 


udd 


1 








A 


uds 


1 





-1 


jr 


sss 


1 


-1 


-3 


7T + 


ud 





+ 1 





7T~ 


du 





-1 





7T° 


{uu + dd) 











K° 


ds 








+1 



Table 1.1: The quark flavor content of some sample hadrons. The properties of the hadrons 
are determined by the properties of the quarks from which they are composed. 




Figure 1.1: The baryon octet (left) and decuplet (right). The N° is the neutron and the 
is the proton. 
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Figure 1.2: The quark model view of hadrons. Left: a baryon consisting of three quarks 
joined by flux-tubes of glue. Right: a meson consisting of a quark-antiquark pair joined by 
a flux-tube of glue. 

The quark model correctly described the baryon octet and baryon decuplet. It also 
predicted that the light meson octet and singlet would mix, forming a meson nonet. This 
corrected the Eightfold Way which had treated a newly discovered meson, the Tj , as a 
singlet with no relation to the existing meson octet. 

Observed states such as the A ++ (uuu) produced from 

7T+ + P + -> A ++ , 

and the Q~(sss) appeared to violate the Pauli exclusion principle, because all three 
fermions were in a symmetric spin-flavor-spatial wavefunction. This led Greenberg to 
postulate the existence of a new quantum number in 1964 [Gre64j which could take on one 
of three values. This additional quantum number later evolved into the concept of 'color,' 
the charge associated with the force holding the quarks together. Thus, every flavor of 
quark comes in three 'colors': 'red', 'green', or 'blue.' The color hypothesis holds that all 
hadrons are colorless states (combinations of all three colors in the case of baryons, or a 
color-anticolor combination in the case of mesons). A profound and disturbing prediction is 
that the fundamental building blocks themselves, the quarks, could never be observed in 
isolation. This phenomenon is known as 'color confinement.' 

Experimental verification for quarks came from the SLAC (Stanford Linear Accelerator 
Center) experiments performed under "Pief" Panofsky |CG89] in the late 1960s using a 
beam of electrons directed onto a hydrogen target. The goal of the experiments was to 
repeat the Rutherford experiment, but at energies high enough (up to about 18 GeV) to 
probe the structure of the proton by deep inelastic scattering. They found that the proton 
did indeed contain concentrations of charge which were 'point-like' in comparison to its 
spatial extent. High energy electrons were unable to knock isolated quarks out of hadrons 
as expected from the color confinement hypothesis. On the other hand, it appeared that 
quarks behaved as free particles within the hadrons. This vanishing of the color force at 
high scattering energies (short distance resolution) is known as asymptotic freedom. 

In 1974, the J/^S, a meson containing a new heavier flavor of quark dubbed the 'charm' 
was discovered independently at Brookhaven and SLAC. Currently, a total of six quark 
flavors has been identified. 
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Other experiments measured the ratio of cross sections 

a(e + e~ — > hadrons) 
a(e+e~ — > 

as a function of center of mass energy. At lower energies, only uu and dd quark pairs can 
be created 1 "!. However, as the center of mass energy increases, it reaches the threshold 
where the more massive ss pairs can be produced, then cc pairs and so on. Because each 
quark comes in one of three colors, there is an overall factor of three in the reaction rate. 
These experiments provided an impressive confirmation of the existence of exactly three 
colors |PRSZ99j . 

~u u~ 




Figure 1.3: The quark model view of the reaction ir~ +p + — > K° + A . Time increases from 
left to right. The up quark and antiquark annihilate, and a new strange quark- ant iquark 
pair appears. 



1.6 The Standard Model 



Currently, the interactions among objects can be understood in terms of three basic forces: 
the force of gravity, the electroweak force, and the strong nuclear force. 

The electroweak and strong interactions are treated in a theoretical framework known as 
the Standard Model, which successfully predicts virtually all observed phenomena in 
particle physics. In the Standard Model, matter consists of quarks and leptons (and 
associated antiparticles), which are described by quantum field theories possessing a local 
gauge symmetry. This local gauge symmetry gives rise to the interactions among the 
quarks and leptons mediated by 'gauge bosons.' The electroweak gauge bosons are the 
massless_7 (photon), and the massive W ± and Z° which acquire their mass (it is 
believed 11 !) from the Higgs boson H via the 'Higgs mechanism.' The strong gauge bosons 
are the massless gluons (see Table 11.21) . 

It is widely recognized that this model is incomplete (most conspicuously, it fails to cleanly 
integrate gravity with the other forces), and there is active research into extensions and 

10 The dynamics of the strong interaction is independent of quark flavor, but the kinematics does in general 
depend on the masses of the different quark flavors. 

11 The Higgs boson has not yet been observed experimentally. 
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revisions. For this work, we will focus on extracting predictions from quantum 
chromodynamics, that component of the Standard Model which defines the dynamics of 
the quarks and which gives rise to the strong nuclear force. 



Generation 


Leptons 


Quarks 


I 


e, v e 


d, u 


II 




s, c 


III 


T, V T 


b,t 



Gauge Bosons 

Table 1.2: The standard model of particle interactions. This model describes all matter 
as made up of quarks and leptons, held together by the interactions of the force-carrying 
mediators. The Higgs boson is the only particle of the Standard Model which has not been 
observed, and is not shown. 



1.7 Quantum chromodynamics 
1.7.1 Quantum field theory 

Quantum field theory (QFT) is the extension of quantum mechanics to relativistic systems, 
allowing for a unified description of matter and radiation fields. The first foundations of 
QFT were laid by Dirac in his 1927 paper [Dir27j which treated the emission and 
absorption of electromagnetic radiation by atoms. 

A quantum field theory is based on a Hilbert space of possible physical states of the 
system. The fundamental degrees of freedom are fields of operators which act on the 
Hilbert space and which are defined over a space-time manifold. 

Natural units are used in which K — c — 1. In Minkowski (flat) space-time, points are 
labeled by ( OC • OC . OC • oc 3 ) = (t, x) with respect to some basis with associated metric 
•q^v = diag(l, —1,-1,-1). Quantities of the forrr0 

V^-y^-v") = (x f "-y")(x fl -y fl ) = (x -y ) 2 -(x-y}-(x-y) 

are invariant under the Poincare group of translations, rotations, and relativistic boosts. 
The fields in QFT may have indices which label different components, and usually these 
components are required to transform irreducibly under the Poincare group [Ram90j . This 
means that an arbitrary Poincare transformation R (e.g. a rotation) transforms the field 
$ a as: 

$ a {x) - & a (x) = ^ b (R- 1 x)D ba (R), 



12 Summation over repeated indices is implied unless noted otherwise. 
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where D(R) is an irreduciblqlj representation matrix representing the effect of R on the 
components of the field. In general, there are many different irreducible representations 
allowing for the definition of a different type of field for each. 

The Poincare group has ten generators. Two operators which commute with all of the 
generators of the Poincare group are called Casimir operators and consist of the mass and 
the relativistic spin. Thus, we label the different types of fields by their mass and spin. 

In the Heisenberg operator picture, we single out the time direction and treat the degrees 
of freedom at a given time as operators with eigenvectors and associated eigenvalues: 

<S>(x;t)\4>;t) = 4>(x;t)\<j>;t). 

For fixed t, all of the spatial points x are space-like separated, which allows us to define 
\4>;t), the simultaneous eigenstate for all of the operators at different spatial locations at 
time t. The eigenvectors represent the possible states of the field at time t and form a 
complete (Hilbert) space: 

(0;t|0';t> = 5(0-0') 
ci0|0;t)(0;t| = 1 

with the appropriate definition of the inner product and integration measure. A general 
state of the system at time t is specified by the state- vector t), and the probability 
amplitudJ^l of the system to be in state \ip'; t') if it was known that it was in the state 
\ip;t) at t < t' is given by the inner product 

A field at a spatial location x evolves in time via the Heisenberg time-evolution equation: 

$(£;£') = e ^ (t '- t) $(x;t)e-^ (t '-* ) 

where H is the Hermitian Hamiltonian operator, the Poincare generator of temporal 
translations. The Hamiltonian operator governs the dynamics of the theory, and its 
spectrum consists of the steady states of the theory, including all stable single- and 
multi-particle energy states. In order for the particle content of the field theory to be 
well-defined, the Hamiltonian must be bounded from below. Subtracting a suitable 
constant from H, we may define the 'vacuum' state \ the time-independent state of 
lowest energy: 

H\Q) = 



13 A representation of a group is a set of matrices {D} which satisfies D(R)D{R') = D(RR') for all R 
and R' in the group. The representation is irreducible if there is no invariant subspace which only mixes 
with itself under the group operations, i.e. if the representation matrices cannot all be simultaneously block 
diagonalized. 

14 In quantum theory, the probability of a particular observed outcome is given by the absolute square of 
the sum of the probability amplitudes for all of the ways that outcome could occur. The fields may not be 
directly observable, but we may always speak of the probability amplitude for a particular configuration's 
contribution to a process. 
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By the action of a suitable operator, any state may be excited from the vacuum: 

\<p ; t) = o^t)\n) 

Because all operators of interest can be expressed as analytic functions of the fundamental 
degrees of freedom of the theory, all of the information about the properties and dynamics 
of the theory is contained in the so called Green's functions, the vacuum expectation values 
of products of fields, such as: 

(n\T[^£ b -,t b )$(s*,t a )]\n), 

where T[- ■ •] is the right-to-left time ordering operator which arranges the fields in time 
order from earliest on the right to latest on the left. 

Inspired by the work of Dirac, Richard Feynman developed the path integral approach to 
quantum mechanics in 1948 |Fey4 8] , which was quickly extended to quantum field theory. 
His work made possible the evaluation of vacuum expectation values (VEVs) by use of the 
Feynman functional integral [PS95j : 

(il\Tl®(x b ;t b )$(x a ;t a )\ \U) = hm hm 



; T^oo(i+i e ) f V(f> exp {z5[0]} 
lim lim — — - — — -- — — rZ\J] 

e-+0T^cx>(l+ie) Z[0] SJ(x b ] t b ) SJ[X a ] t a ) 



J=0 



where 

Z[J)= I Vcj) exp {iS{<p] + %J<j>) 



is a generating functional of the fields, and 



r-T/2 

S[4>] = I dt I d 3 x£(0,<90) 



T/2 JV 

is the action functional, the space-time integral of the Lagrangian density £ which 
determines the dynamics of the theory!^! 

In practice, the infinite-dimensional functional integral over field configurations is 
ill-defined. In order to calculate quantities in quantum field theory, we work in a finite box 
and introduce a regulator which makes the integrals convergent. The regulator is then 
removed by using the method of renormalization, which allows the 'bare' parameters in the 



15 When the magnitude of the action is large (in units of h), then the theory becomes classical and the 
expectation values are dominated by values of the field around which the phase remains stationary. This 
principle of stationary phase (often referred to as the 'principle of least action') gives the Euler-Lagrange 
equations which govern all of classical mechanics: 
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original Lagrangian (such as mass, coupling, and field normalization) to vary as functions 
of the regulator parameter. These bare parameters are not observable, and can be used to 
absorb the divergences encountered in the theory. It is not always possible to do this; in 
non-renormalizable theories the divergences cannot all be absorbed into the small number 
of original bare parameters. 

1.7.2 Reference frame covariance 

The concept of reference frame covariance [Mor83, ISch85] was developed in Einstein's 
Special and General Theories of Relativity and states simply that coordinate systems at 
different space-time points may have different orientations. A physical quantity represented 
by a vector v can always be written as a linear combination of some basis vectors {r}( a )} 
with coefficients v a : 

v = v a rj(a) 

where summation over repeated indices is implied. The basis vectors are abstract entities 
which define some reference frame, and the components are simply numbers which express 
the orientation of the vector within that reference frame. The vector v is an abstract entity 
and is independent of the choice of reference frame. A different reference frame corresponds 
to a different set of basis vectors {v{ a )} anc ^ coefficients {v' a }, but the vector v remains 
unchanged: 

v = v'% a) 

General relativity postulates that a reference frame can only be defined locally. If we want 
to make meaningful comparisons of the components of a vector field {v a (x)} at two 
different points, we need a way to specify the relative orientations of the basis vectors at 
those points. Manifolds are locally flat, which means that we can define vectors connecting 
the points in a small neighborhood of x using a basis {e^)}. Because the neighborhood is 
flat, we can always write the basis vectors at x + ee^) as linear combinations of the basis 
vectors at x: 

fj( a )(x + ee( M )) = f)( a )(x) + igefj{b){x)A^ a {x) + 0(e 2 ) 

which gives the correct behavior as e — > 0. We have pulled out a factor of ig for later 
convenience and have not yet specified the form of the connection A^ix). We now have a 
way of expressing the components of a vector at x + ee^-j on the same basis as the 
components of a vector at x: 

v{x + ee^)) = v a (x + ee^))f)( a) (x + ee^)) 

= v a (x + ee {tl) )(7] {a) (x) + igerj^^A^x) + 0(e 2 )) 

= [{Sab + igtA^ ab (x))v b (x + ee (M) )] r)(b){x) + 0(t 2 ) 

= [U ab {x,x + ee {fl) )v b (x + e (/i) )] fj {a) {x), 

where we have introduced the infinitesimal linear parallel transporter 



U{x, x + ee ip) ) = 1 + igeA^x) + 0(e 2 ) 
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which expresses the components of a vector at x + with respect to the basis at x. 
Alternatively, we may say that the parallel transporter moves a vector from x + ee^ to x 
while keeping it (locally) parallel to its original orientation. 

U(x, x + ee( M )) : (x) <— (x + ee^)). 

We may build up a finite parallel transporter along a directed curve C yx from x to y by 
repeated application of the above. Break C yx up into an N + 1 point mesh {z , Z\, ■ • ■ , zn} 
where zq = x, z^ = y, and the distance between Zk and ^ — 1 is e. Letting dzk = Zk — Zk-i 
we may define the left-to-right path ordered exponential by |Rot97j : 

U(C yx ) = Vexplig I dz" A^z 

I J Cy x 

= lim(l + igdz^A^Zo))^ + igdz^A^zi)) ■■•(! + igdz^A^z^x)) 

We may also define the covariant derivative _D M , which takes into account both the spatial 
change and basis change of v a (x) under an infinitesimal displacement e in the fi th direction 

D v a (x) = lim Uab ^' X + ei ^ V ^ X + e JM} Z V<1 ^ 
- lim ^ ab + igeA vab{%))v h {x + ee (/Lt )) - v a (x) 



(o~abd^ + igA IMjb {x))v b {x) 



The significance of this approach is two-fold. First, the quantity (D^v a (x))fj^ a ) is a proper 
abstract vector; the components transform 'covariantly' under operators on the vector 
space. Second, the connection A(x) encapsulates any non-trivial topological properties of 
the space of different coordinate systems on which the vector is defined. For example, the 
surface of a sphere is locally flat, but has a path-dependent parallel transporter U(C) as 
can be seen in figure 11.41 




Figure 1.4: An example of parallel transport around a non-trivial manifold. The blue arrow 
ends up in a different orientation even though it made only locally parallel moves. 



19 



1.7.3 Quantum electrodynamics 

Hermann Weyl had attempted to apply the concept of local gauge symmetry to 
electrodynamics in 1919. He suggested that one could choose a different scale, or gauge, at 
each point in space, but this idea was quickly proven incorrect (e.g. the atom had an 
observable characteristic scale). He succeeded, however, in the late 1920's along with Fock 
and London to identify the the 'coordinate system' of relevance in electrodynamics. An 
electron's quantum mechanical wavefunction ip(x) can be rotated by a global phase 
without changing the dynamics of the theory. They postulated that the phase of an 
electron's wavefunction was a local quantity and required that comparisons of the phase at 
two different points x and y required the use of the parallel transporter 



U (C yx ) =Vexp\ig f dz» A^z) I 



where the gauge field A^x) is a real- valued field. 

Dirac had already found a quantum field theory describing a free electron when he sought 
to write a relativistic form of Schrodinger equation (which is first order in the time 
derivative) by factoring the relativistic energy-momentum expression |Gri87j : 

GAV ~ m2 ) = (YPv + m)(YPu - m) 
where the four Dirac matrices 7 M must satisfy 

{ff> r} = + 7 V = 2gi*> (1.1) 

The Dirac equation of motion (ry M <9 M — m)^ = comes from the classical limit of the 
quantum field theory described by the Dirac Lagrangian [PS95J: 

£ D irac = ^(h^d^ -m)lp 

where ip(x) is the electron field and ip(x) = ^(x)^ . In order to have a real action (leading 
to a Hermitian Hamiltonian) , we must have (7°7 M )^ = (7°7 M ) and 7 "'' = 7 , or equivalently: 

7 °t = 7 ° ; yt = _y ( L2 ) 



Replacing the regular derivative in ^Dirac whh a covariant derivative gives rise to an 
electromagnetic interaction term: 

il>(x)(i^d f j, — m)ip(x) — > ip(x)(i'y fM D f j i — rn)ij)(x) 

= ip(x)(ij fJ 'd^ — m)ip — gA fl (x) , ip(x)'j IJ 'ip(x) 

-^Dirac ^Int 

Where we see that the gauge field A^x) is to be identified with the electromagnetic 
potential, ipj^ip is the electron probability 'current' (which transforms as a Poincare 
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4- vector), and g = — \e\ is the electron charge, the coupling constant for the 
electromagnetic interaction between the electron field if) and electromagnetic field A. 

This theory possesses an important symmetry: local gauge invariance. The Lagrangian is 
unchanged under local gauge transformations (picking a new phase at each space-time 
point, and updating the connection accordingly): 

V>(x) -> e ta{x) ^(x) 
^0) -> ^(x)e" iQ(x) 

A M (x) -> A^x) - -d^afa). 

Quantum field theories which are based on a local gauge symmetry are known as gauge 
theories. In 1971, Gerard 'tHooft [tH71J showed that all gauge theories were renormalizable. 
This was a critical success for gauge theory, because it means that calculations are 
guaranteed to give finite results for physical processes as the regulator is removed. 

In analogy with General Relativity and figure 11.4} the field strength tensor is associated 
with the curvature of the gauge field, the amount the phase changes around an 
infinitesimal square loop pch8 5, MM94) in the /x — v plane: 

F^ = -[D^,D U ] 
W 



Considering first 



D^Dv = —(0^ + igAp{x)) (d u + igA v {x)) 

= —d^ + duA v (x) + A v (x)duAJx)d v + ig 2 A ll A v 
W 



we see that 



F,j, v = —[Dim, D u ] = d^A u -d u A^ + ig[A^A u ] (1.3) 
W 

= D^A U — D u A fjL 

The gauge group of electrodynamics is the group of all unimodular complex phases U(l). 
Because U(l) is Abelian (commutative), the commutator in 11.31 vanishes and we are left 
with the familiar expression 

Fp,v = d^A u — dyA^ 



The kinetic term in quantum electrodynamics gives Maxwell's equations in the classical 
limit, and is obtained by considering the simplest gauge- invariant quantity involving the 
field strength tensor: 

■^Maxwell T D^yF^ 1 . 
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Putting it all together gives the Lagrangian describing quantum electrodynamics (QED), 
one of the most successful theories in the history of physics: 



£qED — ^Dirac + Ant + ^Maxwell 

1 

4 

Dp = dp + igAp 



Fp U = DpA v — DyAp — dpA v — d u A^ 



QED is formulated on a flat space-time manifold, but may have curvature in 'phase space' 
described by A^(x). If the phase space is flat, then the parallel transporter will be path 
independent, meaning that A is a total derivative of some function. In such a case, we can 
perform a local gauge transformation which eliminates the electromagnetic field: A — ► 0. 
Thus we see that if we treat the electron phase as an internal degree of freedom, it can give 
rise to electromagnetic interactions as a result of curvature in this internal coordinate space. 

The gauge field can be quantized to give the mediating boson of the gauge force. U(l) has 
only one parameter A^(x), and thus there is only one gauge boson of the quantized theory, 
the photon. 



1.7.4 From pions to the Standard Model 

In 1936, Breit |BCP36] discovered that the strong nuclear force was blind to electric 
charge. That is, the strong nuclear force between two protons was the same as the strong 
nuclear force between a proton and neutron, or between two neutrons. Heisenberg 
suggested that the proton and neutron were two states in the 'isospin' doublet: 




Yang and Mills promoted isospin to a local gauge symmetry by introducing a connection 
based on the SU(2) gauge group, the set of all complex 2x2 matrices with unit 
determinant. This allowed for the identity of the proton and neutron to vary as a function 
of position. The Yang-Mills SU(2) gauge field was shown to describe the pion in the same 
way that the U(l) field describes the photon. In the case of SU(2), there are three 
generators yielding three gauge bosons, the 7r~, 7r°, and n + . 

Physicists then searched for the 'correct' gauge theories describing the interactions of the 
standard model. Work on the electroweak Lagrangian describing the interactions of leptons 
and quarks via the W ± , Z°, and 7 gauge bosons was started in 1961 by Glashow [Gla61j . 
completed in 1967 by Weinberg |Wei67j and Salam [Sal68j . and is based on the 
SU(2) x U(l) gauge group. 



22 



1.7.5 Quantum Chromodynamics 

Quantum chromodynamics (QCD) is the gauge theory describing the interactions of quarks 
and is the primary focus of this dissertation. Fritzsch and Gell-Mann |FGM721 IFGML73| 
identified the quark color as the fundamental charge associated with QCD. Quark fields 
come in six flavors, possess four spin components (which transform irreducibly under 
Poincare transformations), and have three color components. Let the quark field be 
denoted by ipaa( x )y where A is the flavor index, a is the spin index, and a is the color 
index. In this section we will focus on the behavior of the color index, and will suppress 
indices where possible for simplicity of notation. 

The color reference frame (defining the directions of 'red', 'green', and 'blue') is defined 
locally. The infinitesimal linear parallel transporter defining the covariant derivative is now 
an SU(3) matrix instead of a U(l) phase: 

A c 

U ab (x, x + ee (At )) = exj>(igeAl(x)—)ab 

where the A a are the eight 3x3 traceless Hermitian Gell-Mann matrices which generate 
the SU (3) group. SU (3) is the non-Abelian (non-commutative) group of 3 x 3 complex 
unitary matrices with unit determinant. 

We get the covariant derivative = + igA^ by defining 

This gives us eight real gauge fields A a {x) which transform as vectors under Poincare 
transformations. They represent eight spin-1 gluons, the gauge bosons which mediate the 
QCD interactions among quarks. The spin-1/2 quarks are represented by Grassmann 
(anticommuting) fields which transform as Dirac spinors under Poincare transformations. 
In a hypothetical free field theory where A = (all the links are unity), the quantum 
operator \l/ creates a quark state and annihilates an anti-quark state while \& annihilates a 
quark and creates an anti-quark. In the general interacting theory (A^O), we will still 
refer to \I/ as a quark source and \1/ as a quark sink. 

The gluon field strength tensor is given again by the curvature of the gauge field, this time 
with a non- vanishing gauge field commutator [A^, A v \. 

F^ = -[D^,D U ] = D^-D^ 

= d^A v - d v A^ + ig[A^A v ], 

The new commutator term can be expressed in terms of the gluon fields by using the fact 
that the Gell-Mann matrices {A a } obey the SU(3) Lie Algebra 



[A tt , Afe] — i2f a b C X c , 
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where the f abc are the completely antisymmetric SU(3) structure constants defining the 
group. 

ig[A„A v } = |^4[A a ,A 6 ] 

= — ^fabcA. a ^A h v X C = —-fabcA^A^X 11 

Therefore 

F = -F a \ a 

FZv = d,Al-d u Al-gf abc A b A c 

Using the fact that Tr(A a A fe ) = 26 ab , we may put these gauge-invariant pieces together to 
define the QCD Lagrangian density, the functional defining the dynamics of the interaction 
of quarks and gluons: 



QCD 


= iPWDp - M)V - ^Tr(F^F^) 




= 8^ + igAp, 




= -A a u \ a 
2 M 




= 2 ^ ' 




= d^A a v -d u Al-gf abc A b A c 


M 


= di&g(m d ,m u ,m s ,m c ,m b ,m t ). 



In the above, M in a 6 x 6 matrix in flavor space, 7 M is a 4 x 4 matrix in spin space, and A' 
is a 3 x 3 matrix in color space. In other spaces, the matrices appear as the identity (e.g. 
7 M acts as the identity in color space). The trace is taken with respect to color. The 
non-Abelian nature of QCD introduces three- and four gluon interactions through the 
kinetic term — Tr(i 7 ) u ,i 7V " y )/2. These gauge field self-interaction terms are believed to 
reproduce the desirable properties of asymptotic freedom and color confinement. Also, the 
strong nuclear force binding nucleons together in the nucleus is a seen to be a van der 
Waals-type residual interaction among the quarks in different nucleons. 



1.7.6 Calculating in quantum chromo dynamics 

The Lagrangian density £qcd completely determines the dynamics and properties of 
quantum chromodynamics. However, theorists still need to calculate complicated integrals 
of the fornj^l 

i J V^VAf(iJ^,A)exp{iSg,^,A}}, 



The exact integrals we will need will be discussed in the next chapter. 
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where 

Z = J V^V^VA exp {iS$, if), A]}. 

At high energies, asymptotic freedom tells us the the interactions among quarks are weak, 
and can thus be described using perturbation theory. Perturbative QCD treats the 
exponential as a power series expansion in the strong coupling constant g. At high 
energies, or equivalently short length scales (such as those probed in deep inelastic 
scattering experiments), the renormalized coupling constant becomes small and the power 
series is convergent. By calculating the 'tree level' and the first few 'correction' terms for 
various processes, theorists have shown that the predictions of QCD are in striking 
agreement with the results of high-energy experiments [Mut98j . 

At lower energies (e.g. below 2 GeV in e+p scattering, where Bjorken scaling is violated), 
we enter the domain of 'nonperturbative phenomena,' such as color confinement 
(hadronization) . At these energies, the renormalized coupling constant becomes large, and 
perturbation expansions diverge badly. 



1.7.7 The lattice regulator 

In 1974, Kenneth Wilson introduced a regularization method |Wil74bj which made rigorous 
low-energy calculations in QCD possible for the first time. Wilson defined an action on a 
discrete space-time lattice with spacing a. This serves two purposes: first, if the lattice is 
introduced in a finite volume, the number of integration variables becomes finite and the 
integrals become well defined. Second, the lattice regulator introduces an ultraviolet cutoff 
at the scale of the inverse lattice spacing 1/a, thus controlling divergences. The 
introduction of a regulator necessarily breaks symmetries of the theory, and it is hoped 
that these symmetries will be recovered in the continuum limit. For instance, the lattice 
regulator explicitly breaks Poincare invariance. 

Let the elementary lattice vectors which give the minimal displacement in each direction be 
denoted as: 

p, = ae^). 

The quark fields are defined at each site as 

where x denotes the discrete site on the lattice, A is the flavor index, a is the Dirac spin 
index, and a is the color index. Rather than simply replacing the derivatives in the 
continuum action with finite differences, Wilson preserved local gauge invariance by 
working directly with gauge links, the SU(3) parallel transporters defined on the 
elementary links between neighboring sites: 



U^ab(x) = U ab (x,x + p) 
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where a and b are color indices. For simplicity of notation, we will continue to suppress 
indices where possible. SU(3) matrices U and Q will always act on the color degrees of 
freedom, Dirac gamma matrices 7 M will always act on the spin degrees of freedom, and the 
mass matrix M will always be diag(mrf, m u , ■ ■ • ) and will act on the flavor degrees of 
freedom. 

Under a local gauge transformation represented by SU(3) matrices Q(x) at each site, the 
quark fields and gauge links transform as: 

ijj(x) — > tt(x)if)(x) 

ip(x) — > ip(x)QJ(x) 

U^x) -> n(x)U^x)Q^(x + p) 

Ul(x) -> n(x + fi)U^x)tt\x) 

This formulation allowed Wilson to form gauge-invariant expressions such as 

ip(x)U lx (x)ip(x + p) 

and 

tp(x + p,)Ul(x)i/)(x) 

denoted pictorially in Figure 11.51 

U,(x) U\{x) 
ip(x) • < • ip(x + ji) ijj(x) • > m ijj(x + pL) 

Figure 1.5: Two gauge-invariant quark- antiquark operators. The quark and antiquark fields 
at neighboring sites x and x + fi may be combined by use of a gauge link at x. The gauge link 
Ufj,(x) parallel transports a color vector from x + fx to x and the Hermitian conjugate gauge 
link W(x) transports a color vector in the opposite direction. Any gauge-invariant quark- 
antiquark operator can be formed by connecting the fields at any two sites by a suitable 
product of gauge links. 

If we replace D ll il){x) with the symmetric difference 

-V^z) = g- [Up{x)il)(x + fi)~ U\{x - p)ifi(x - ft)) 
we get the fermion part of the action (in Minkowski space-time): 

Sjf = ^a 4 ^(x)(* 7 M -V M -MMo;) 

x 

The exact form of fermion action we will use will be discussed in chapter [2j It differs from 
the above because we will work in Euclidean (not Minkowski) space-time, use an 
anisotropic lattice, absorb the factors of a to work in 'lattice units,' and add the 'Wilson' 
term which fixes a discretization problem (the 'fermion doubling' lattice artifact). 
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Another gauge-invariant quantity involving just the links is the plaquette, the trace (with 
respect to color) of a product of links around an elementary loop on the lattice: 

= Tr(U^x)U u (x + fi)Ul(x + u)U\x)), 

denoted in figure 11.61 

X + V 




X X + jl 

Figure 1.6: The gauge-invariant plaquette operator U^ u (x). 

The plaquette is the trace of the parallel transporter around an elementary loop. It is 
reasonable to expect that the plaquette is related to the field strength tensor, and this can 
be verified by writing U^ v (x) in terms of the A fields using the Campbell-Baker-Hausdorff 
formula: 

exp(v4) exp(I?) = exp(A + B + - [A, B] + (terms higher in A and B)). 

Because we are interested in the continuum limit where a — > 0, we may write 

Uy,{x) a ^>° expligaA^x)) 



Then expanding the exponentials and keeping terms only up to order a 2 gives: 

CV(x) ee U„{x)U v {x + fl)Ul(x + v)U\x) 

= exp(igaA^(x)) exp(igaA u (x + /t)) exp(— igaA^x + v)) exp^—igaA^^x)) 
= exp (tga 2 (VlA u (x) - V s u A^x) + ig [A„(x), A u (x)\) + 0(a 3 )) 

^ exp (iga 2 F )lu (x) + 0(a 3 )) 



where 



1 



V f u A v (x) ee ~{A v {x + p)- A v {x)) a ^ d,A u 



x 



a 

is the forward lattice difference when applied to gauge fields, and we have evaluated all of 
the commutator terms at x, for example: 

-^[A^Mx + ft] = - 9 ^r\A,(x),A u (x)]+0(a 3 ) 
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Thus 

ReTr(CV) = I(Tr(^) + Tr(f/tj) 

= ^(Tr(fV) + Tr(f^)) 

= Tr(l) -^V^Tr(F^F^) +0(a 5 ) (no sum) 

because F u/M = —F^ v (also because Tr(F Ml/ ) = for 577(3)). In the gauge action, we will 
write iV = Tr(l). 

The pure gauge (also called Yang-Mills) part of the lattice action (in Minkowski 
space-time) is thus given by [MM94] : 

S% = -|^5]ReTr(l-^)), /3 = ^, 

= -2EE a4 ^(vn+o(a 5 ). 

For SU(3), P = p-. In the case of Abelian gauge groups, the expression for /3 is divided by 
two in order to reproduce the continuum action to leading order in a [MM94J. 

1.7.8 Renormalization and integration 

Wilson's lattice regulator breaks Poincare invariance, but it maintains local gauge 
invar iance which is crucial to ensuring that the regularized theory is renormalizable. In the 
continuum limit (lima — ► 0), it is easy to see that Poincare invariance is restored. Any 
quantity calculated using the lattice regulator will equal the continuum quantity plus an 
infinite number of 0(a) terms. There is no way to guarantee that the total effect of these 
terms will vanish as a — > (a series of terms may combine to give an 0(1 /a) contribution 
which will remain in the continuum limit). Thus it is critical to make a theory which is 
renormalizable at finite a. 

Another advantage is that the Wilson action is a functional of the gauge links U^(x), which 
are compact gauge fields. This means that the parameter space we must integrate over is 
closed. In contrast, the underlying gauge fields A a (x) may take any value in the open set 
(— oo, oo). There are an infinite number of values of A^x), related by gauge 
transformations, that correspond to a single value of the gauge link U^(x). This gauge 
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freedom implies that gauge fixing must be used in integrals over the A fields to avoid 
divergences. No gauge fixing is required when integrating over the gauge links 
directly |Wil74bj . 

The coupling constant g is renormalized by tuning (3 = Q/g 2 and calculating some physical 
quantity to get a. In this way, we may send a to zero by increasing (3. When all ratios 
become independent of a, then we have reached the scaling regime. Similarly, we may 
renormalize the quark masses by tuning the bare quark masses while measuring the pion 
mass m v and watching the ratios of different calculated masses. It is hoped that these 
parameters will 'flow' to a critical point where the theory will reproduce the predictions of 
continuum QCD. 

The appearance of the inverse power of the coupling constant in j3 = Q/g 2 allowed Wilson 
to calculate some, but not all, non-perturbative quantities using the so called strong 
coupling expansion. For instance, he showed that the static quark-antiquark potential 
increases linearly with separation in the Yang-Mills (pure gauge) SU(3) theory, which is 
evidence that QCD may be confining in the continuum limit. 

The lattice regulator is also perfectly suited for calculation on a computer. The framework 
of calculating QCD quantities using a lattice regulator on a computer is called lattice QCD. 
Starting from the work of Michael Creutz [CJR79] and others in 1979, lattice gauge 
theorists have applied the techniques of Monte Carlo integration discussed in chapter [3] to 
calculate increasingly ambitions quantities. This work is a contribution to one such effort: 
quantum operator design for the low-lying baryon spectrum. 



1.8 Unresolved mysteries in hadronic physics 

The high energy predictions of QCD obtained from perturbation theory do a spectacular 
job of describing the measurements of high-energy physics experiments. On the other 
hand, the low-energy behavior of QCD has not been thoroughly tested. This is due to the 
fact, mentioned above, that such calculations are nonperturbative with respect to the 
coupling constant of the theory. Fortunately, theoretical methods and computational power 
have now evolved to the point where the extraction of low-energy predictions from QCD 
has become feasible. 

Predictions from lattice QCD will come at an opportune time; there are several 'mysteries' 
concerning the hadron spectrum which the constituent quark model, along with other QCD 
inspired models, have failed to resolve in a definitive manner. 
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1.8.1 Dynamics in the constituent quark model 

The constituent quark model is able to predict the masses of many observed low-energy 
hadron resonances. It does this by treating the hadron states as either a quark-antiquark 
state (for mesons) or as a three-quark state (for baryons). The quark interactions lead to 
predictions of the total energy and quantum numbers of the state. Different quark models 
specify different interaction mechanism, but most consist of non-relativistic dynamics of 
the form (adapted from |H + 02j ): 

1. A confining interaction, which is generally spin- independent. 

2. A spin-dependent interaction, modeled after the effects of gluon exchange in QCD. 
For example, in the S-wave states, there is a spin-spin hyperfine interaction of the 
form 

H HF = -gMj2(° X aU° X a)j 

i>j 

where M is a constant with units of energy. Spin-orbit interactions, although allowed, 
seem to be small. This may be due to relativistic QCD-inspired corrections such as a 
Thomas precession term [CI86J. 

3. A strange quark mass somewhat larger that the up and down quark masses, in order 
to split the SU (3) flavor multiplets. 

4. In the case of isoscalar mesons, an interaction for mixing qq configurations of different 
flavors (e.g. uu <-> dd <-> ss), in a manner which is generally chosen to be flavor 
independent. 

1.8.2 Missing baryon resonances 

The constituent quark model predicts many more baryon resonances than have been 
observed experimentally. Because most of the observed baryon resonances have been 
produced from electron beams (electroproduction) or photon beams (photoproduction), it 
is possible that these 'missing resonances' are created in channels which have not yet been 
sufficiently examined. Experiments are currently underway using different types of beams 
(such as linearly polarized photons) to excite resonances through different processes. 

An alternative explanation suggests that the missing resonances are due to an excess of 
degrees of freedom in the quark model. One proposed fix is to reduce the degrees of 
freedom by 'locking' two of the quarks of a baron together in a di-quark configuration. In 
this picture, the higher resonances would be various single-quark excitations of the 
three-quark bound state. 

Direct QCD calculations are needed to determine the properties of these missing 
resonances. Lattice QCD is in a position to identify the relevant degrees of freedom in a 



30 



Name 


L 


s 


P 


Oscillator Band 


iV(939) 





1/2 


+ 





iV(1440) 





1/2 


+ 


2 


iV(1535) 


1 


1/2 




1 


iV(1650) 


1 


3/2 




1 


JV(1710) 





1/2 


+ 


2 



Table 1.3: The lowest lying spin- 1/2 nucleons. The numbers in parenthesis are the masses 
of the resonances in MeV, L is the orbital angular momentum of the quarks, S is the total 
spin of the quarks, and P is the parity of the state. The quark model treats the nucleon 
resonances as increasing excitations in an oscillator potential. The iV(1440), or Roper, is 
measured to be lower than predicted by the quark model. 

baryon. Simulations can probe the spatial structure of the fields which compose the baryon 
state thus leading to a better understanding. 

1.8.3 The Roper resonance 

The simplest confining interaction in the quark model is the linear oscillator potential. If 
that potential is used, then the baryon spectrum can be divided into oscillator bands, 
which correspond to increasing excitations of the spatial wavefunctions with alternating 
parity. Within each band, the levels split according to spin-spin and spin-orbit interactions. 
Consider the five lowest known spin- 1/2 nucleons (J=l/2, 1=1/2) and their corresponding 
levels in the oscillator approximation |H + 02] shown in table 11.8.31 

The iV(1440), known as the 'Roper' resonance, is the first excited even-parity nucleon 
state, but it is below the lowest odd-parity spin- 1/2 nucleon state, the iV(1535). This is 
surprising |MCD + 05] , because most quark models treat the Roper as a radially-excited 
'breathing mode' of the iV(939), which belongs to the second oscillator band having two 
nodes in the spatial wavefunction. On the other hand, the odd-parity state iV(1535) is 
modeled as an orbitally excited state in which the quarks have orbital angular momentum 
L — 1 with only one node in the spatial wavefunction. 

Alternative models have proposed that the Roper resonance is a qqqg hybrid state with 
excited glue or that it is a qqqqq five-quark molecular state. The mystery is complicated by 
the facts that the Roper has a wide width and has only been detected from phase shift 
analysis, never directly as an enhancement in the cross section. 
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1.8.4 The A(1405) 

The A(1405) having I(J P ) = 0(| ) is another puzzle in baryon spectroscopy |NNMS04] . In 
the quark model, the A series of resonances are described as different excitations of a 
three-quark uds flavor singlet system. Because the A contains a strange quark, the first 
odd-parity state would be expected to be above the first odd-parity nucleon, which contains 
only the lighter up and down quarks. However we see from table [TA31 that the A(1405) is 
below the iV(1535), the first negative-parity nucleon. It has been proposed that the A(1405) 
is a linear combination the three-quark uds state with two baryon-meson bound states NK 
and 7rS, where the S baryons have even parity, strangeness=-l, and I{J P ) = 1(| + )- 

Because the molecular states will mix with the three-quark state, lattice calculations are 
needed to determine the amount qqqqq contribution vs the amount of qqq contribution in 
the A(1405). Understanding molecular hadronic states will lead to a better understanding 
of hadron structure and hadron reactions. 



1.8.5 Exotic hadron states 

The quark model treats hadron resonances as excitations of quarks in a confining potential 
with inter-quark interactions. It does not have anything to say about excitations of the 
gluon field. Laboratories are dedicating significant resources to the exploration of these 
'hybrid' states of quarks bound by excited glue. QCD predicts that the gluon field can 
carry quantum numbers such as angular momentum, resulting in the possibility of exotic 
states, with quantum numbers inexplicable within the quark model. If such a state was 
found, it would be a significant low-energy validation of QCD. Two exotic J PC = 1 h 
hybrid meson candidates have been tentatively identified at 1.6 GeV and 1.4 GeV in the pn 
and r/ir channels respectively by the E852 collaboration at BNL [A + 98l lC + 99] . 

Lattice QCD theorists are assisting in this effort by calculating the properties of exotic 
mesons and baryons. Understanding hybrid states will lead to a better understanding of 
when and how the gluonic degrees of freedom play a role in the baryon spectrum. 

1.8.6 Current experimental efforts 

The spectrum is a fundamental quantity of QCD, yet remains incompletely understood by 
experimentalists and theorists alike. There are a number of high-precision baryon 
spectroscopy initiatives at laboratories including Jefferson Laboratory (JLab), BNL, 
Mainz, Graal, and BES. 

CLEO-c at Cornell is committed to exploring heavy quark systems |Col01j . Hall B at JLab 
is dedicated to mapping out the low-lying baryon spectrum, and the planned Hall D at 
JLab will search for exotic hadrons [Man05] . Each one of these efforts will rely on lattice 



32 



QCD calculations for valuable theoretical input. Such calculations will tell the 
experimentalists the quantum numbers of the states they may find, and will help them 
decide on the optimal energy ranges to explore. 



1.9 Goal of this work 

The need has never been greater for a systematic exploration of the hadron spectrum using 
the best theory we have, quantum chromodynamics. Unfortunately, low-energy features of 
QCD cannot be treated by perturbation expansions in the coupling constant. As the 
world's experimental dataset grows, so does the number of resonances which cannot be 
adequately described by QCD-inspired phenomenological models. The lattice approach 
allows us to numerically solve QCD, not just model it, and is currently the best tool for 
low-energy QCD calculations available. 

The Lattice Hadron Physics Collaboration (LHPC) was formed in 2000, with the stated 
goal of using lattice QCD to understand the structure and interactions of hadrons [l + 00j . 
One major goal of the collaboration is to calculate a significant portion of the low- lying 
spectrum of baryon resonances in QCD by developing and using the best possible methods 
in lattice QCD. 

Spectrum calculations serve four main purposes. First, such calculations help answer the 
question: can QCD reproduce the observed hadron spectrum? Second, predictions of the 
energy and quantum numbers of new states will help the experimental community as they 
plan and design the next generation of experiments. Third, spectrum calculations help 
identify the important degrees of freedom in the theory. Specifically, theorists hope to 
identify subsets of configurations which dominate the functional integral. In 
thermodynamics, systems with an unmanageably large number of degrees of freedom can 
be described in terms of a few bulk properties such as pressure and specific heat. Similarly, 
lattice QCD may find a way to describe systems in quantum field theory in terms of a few 
dominant subsets of configurations. Such a description would facilitate the construction of 
simpler models and calculation schemes which would serve as an analytical complement to 
a fundamentally numerical approach. Among other things, such models may help us better 
understand the process of hadron formation in the early universe. Fourth, just as hydrogen 
spectroscopy led to the rapid development of quantum mechanics, so too may hadron 
spectroscopy lead to new theoretical breakthroughs. 

This work presents the operator design, tuning, and pruning techniques developed during 
the course of my research. These methods have been used to achieve unprecedented signals 
for spectral states in the nucleon channel, and are readily extended to the other baryon 
channels (A, A, E, S, • • - ). This work describes both the methodology used in the operator 
design and the computational steps involved in tuning and pruning the operators. Specific 
operators are presented which can be used to extract the low-lying nucleon spectrum in 
future high-statistics, low-quark mass calculations. Preliminary heavy-quark quenched 



33 



spectrum results are given which, while unsuitable for comparison with experiment, do 
show the signal quality achieved by the operators. It is expected that these operators will 
play a valuable role in the LHPC's ongoing high-precision baryon spectrum calculations. 



1.10 Organization of this dissertation 

The formalism of extracting the spectrum of a quantum field theory is reviewed in chapter 
[2J We describe how QCD is formulated on a computer, how the Monte Carlo method is 
used to calculate expectation values of correlation functions of quantum operators, and 
what these functions can tell us about the spectrum of QCD. We discuss the role of 
Euclidean space-time correlation functions, a method for reliably extracting multiple 
excited states, and the importance of the lattice symmetry group in identifying the spin of 
the created states. 

In chapter [3], we introduce the lattice-regulated formulation of QCD, and the review the 
Monte Carlo method for evaluating functional integrals. We also discuss our error analysis 
methodology which uses bootstrap and jackknife resampling. 

In chapter HI we discuss construction of nucleon operators from the building blocks we have 
at our disposal on the lattice. We introduce smeared quark- and gauge-fields which provide 
a closer link between our operators and phenomenological models. Covariant displacement 
is discussed as a natural way to build up radial structure. We treat in detail the 
classification of our quantum operators based on their behavior under groups of symmetry 
transformations, and how this behavior relates to the continuum spin values of the states 
created by our operators. 

In chapter [5j we discuss the combination of our operators into correlation matrices. We 
integrate out the quark degrees of freedom and reduce the problem to that of evaluating 
three-quark propagators. We then provide a detailed walk-through of the computational 
steps performed on the Carnegie Mellon University Medium Energy Physics computer 
cluster which allowed us to tune and select the best quantum operators for extracting the 
low-lying QCD spectrum. 

The tuning of the operators is discussed in chapter [6l where stout-link gauge-field smearing 
is shown to reduce the noise of extended baryon operators, while Laplacian quark field 
smearing is shown to reduce the contamination in of our baryon operators due to 
high-frequency modes. The systematic approach we developed to tune the smearing 
parameters is presented, and the optimal smearing parameter values are reported. 

The construction of nucleon operators utilizing the symmetries of the lattice discussed in 
chapter 0] leads to a unmanageably large number of operators for our correlation matrices. 
Chapter [7] discusses the systematic approach we developed to prune these operators down 
to an optimal linear combination of sixteen operators. 
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These sixteen operators are then used to extract the low-lying nucleon spectrum on 200 
quenched configurations. In chapter [S] the methodology and analysis code written for this 
task are discussed, and fit results are presented. Although the configurations are quenched, 
in a small volume, and have an unrealistic quark mass, we still get reasonable results for 
the nucleon spectrum. The spectrum extracted in this work is discussed in comparison to 
the experimentally measured nucleon excitation spectrum, and the predictions of the 
relativistic constituent quark model. We conclude by summarizing the main results of this 
work, and by discussing the ongoing efforts of the Lattice Hadron Physics Collaboration. 



Chapter 2 



Calculation overview 



2.1 Spectral states and resonances 

2.1.1 Hilbert spaces 

In quantum field theory, all physical states of a system are represented by rays in a Hilbert 
space Ji. A ray is an equivalence class of vectors differing only by a non-zero scale factor. 
Thus |0) and a\(p) (a^O) represent the same state of the system. 

A Hilbert space is defined as a normed (or metric) space in which all Cauchy sequences are 
convergent. The norm || • || is defined for any vector |0) as 

II 10) IN \fW)- 

A sequence of vectors (|0 n )} is a Cauchy sequence if for every real number e > there is a 
non-negative integer N such that 

II \<t>n) - \<t>m) II < e V m,n> N. 

A normed space is a Hilbert space if all Cauchy sequences in the space converge to an 
element in the space. This property is essential to justify many assumptions we will make 
about the convergence of integrals over elements in this space. 

2.1.2 The Hamiltonian 

Allowed transformations on the system are represented by unitary operators which map 
Hilbert space states to Hilbert space states. Examples of transformations include Poincare 
transformations (space-time translations, rotations, and boosts), reflections, and charge 
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conjugation. Continuous transformations can be written as the exponentiation of 
Hermitian operators, known as generators. For example, spatial translations are generated 
by the momentum operator p: 

f(x)^f'(x) = f(x-a) 

= exp(-a 3 dj)f(x) 

= exp(ia J pj)f(x), pj = idj 

temporal translations are generated by the Hamiltonian, and rotations are generated by 
the angular momentum operator. 

The spectrum of a theory consists of those states \k) in the Hilbert space which are 
eigenstates of the Hermitian generator of temporal translations H, the Hamiltonian 
operator: 

H\k) = E k \k). 

These are the steady-states of the theory which remain stable if allowed to freely propagate. 
Because the Hamiltonian operator is Hermitian, the eigenvalues are real and represent 
the energy of the states \k). States which are eigenvectors of any Hermitian generator 
which commute with the Hamiltonian (such as linear and angular momentum) may also be 
labeled by their associated eigenvalue. Thus, it is possible to speak of a state having a 
particular value of both energy and momentum. 

An additional constraint placed on a quantum field theory is that the Hamiltonian must be 
bounded from below. This means that there is a stable state of minimum energy, the 
vacuum \ fl). By adding a suitable constant to the Hamiltonian (which does not change the 
spectrum), we may require that 

H\Q) = 0. 

The existence of a state of minimum energy is essential in a quantum field theory to ensure 
that the field degrees of freedom in a finite spatial region don't radiate away an infinite 
amount of energy as they drop to lower and lower energy states. 

The stable energy states \k) form a complete basis for the Hilbert space: 

5>><*l = l, 
k 

where the sum is over both discrete and continuous states (with the appropriate measure). 
2.1.3 Resonances 

In a scattering experiment, a resonance is characterized by a dramatic increase in 
cross- section with a corresponding sudden variation in the scattering phase shift. A 
resonance observed in the scattering of hadrons can often be interpreted as an unstable 
particle lasting long enough to be recognized as having a particular set of quantum 
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numbers. Such a state is not an eigenstate of the Hamiltonian, but has a large overlap onto 
a single eigenstate. Note that if the quark mass is taken to be large in the simulation, then 
it may happen that a normally unstable state would be kinematically forbidden to decay 
due to energy conservation. Because its would-be decay products are too heavy, such a 
resonance would appear as a stable state in the spectrum. 

Finite spatial volumes provide a second way to infer the existence and properties of 
resonances. A method of using the lattice approach to probe resonances is presented in 
Refs. |Lus86[ ILus91bl ILus91aj . In this method, the dependence of a two-particle spectrum 
on the spatial box size is used to determine elastic scattering phase shifts, which are in 
turn related to the energies and widths of any resonances in that channel. For this work, 
we will extract the spectrum in a single volume. Future studies will be needed to determine 
the single and multi-particle content of the spectrum we find. 



2.2 Euclidean space- time 

Having the Lagrangian describing QCD, our next step is to find a method to extract the 
spectrum. In ordinary quantum mechanics, we would simply solve the Schrodinger 
equation for the stationary states in time: idtip = Eip. Unfortunately, it is extremely 
difficult to work directly with the QCD Schrodinger equation due to the number of 
variables and the self-interaction terms. It is more practical to extract the energy levels of 
the theory from correlation functions of fields in Euclidean space-time. 

The following discussion will deal with the theory of a single real scalar field 0. Afterward, 
we will discuss subtleties that arise when applying the results to the case of quark fields 
and gauge links. 

The Hilbert space and the original field can be reconstructed from the Green's functions. 
Osterwalder and Schrader |OS73al IOS75] showed that under general conditions, we may 
analytically continue the Minkowski Green's functions to imaginary time. Doing so will 
make the integrals explicitly convergent and will allow us to use the computational tools of 
statistical mechanics to calculate quantities in this system. We will summarize the main 
features of this approach in the following. A more detailed treatment is presented 
in |MM94j . For the rest of this work, we will work in Euclidean space-time, and will denote 
Minkowski space-time quantities with a subscript or superscript M. 

The points on the Euclidean space-time manifold are related to those on the Minkowski 
space-time manifold as: 

— ju 4 — ]\/[ — - ju — J — iV/ — j ' * ' 

d 4 = d 4 = -id™ = -idP M , d 3 = cP = df = -ff M . (2.2) 

With these definitions, the new metric therefore becomes Euclidean: 

5^ u = diag(+l, +1, +1, +1), and there is no distinction between covariant (lower) and 
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contravariant (upper) indices. Summation over repeated indices will be assumed unless 
specified otherwise. 

Making the above replacements, the continuum fermionic action becomes: 

iSp = i [ dx M d 3 x M ip M (x M )(ilMdo +i r f M df -m)ip M (x M ) (2.3) 

dxid 3 xip(x)(— 7°/<9 4 + i^fijdj — m)if){x) (2.4) 
If we define the Euclidean space-time Dirac 7 matrices as 

7 4 = 74 = TXr, Ik = 7" = -il k M (2-5) 

then we have: 

iS% -> - J d 4 x^(x)(^ + m)ij(x) (2.6) 

= -S F . (2.7) 

Using Eqn. 11.11 and Eqn. 11.21 we see that the Euclidean space-time Dirac 7 matrices 
defined in 12.51 satisfy: 

{7m, 7*} = 2^, 7l = 7„ (2.8) 
We will use the standard Dirac-Pauli representation for the 7 matrices: 



where J is the 2x2 identity matrix and the Pauli spin matrices are 

•>-(:;). -(ro 4 ). -(;_ o j. ^ 



The fermion fields ip{x) and if)(x) transform irreducibly under SO '(4) rotations in Euclidean 
space-time, and represent the Minkowski fields analytically continued to imaginary time. 
Unlike in Minkowski space-time, if) and ip must be treated as independent fields: if) 7^ V^74- 
This is required to simultaneously satisfy Euclidean covariance of the fields, the canonical 
anticommutation relations of the fermion fields, and the equality of the Euclidean 
two-point function with the relativistic Feynman propagator continued to imaginary 
times |OS73bllWil74aj . 
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2.3 Euclidean gauge links and gauge action 



The fermion action may now be discretized as before, using gauge links which are now 
defined between neighboring quark field sites on the Euclidean manifold. 



X+jJ, 



Vexp l^ig J dz v M A v (z 
Vexplig / dz v A v (z) 



with the Euclidean gauge field defined by 

A^ = A = ^Aq = iAj^^ Aj = A 3 ' = Aj = A^ 

This gives 



m- 



F*> - -F w - -F n 



— or M — ir Q j , 



pjk _ p. 



<jk _ t?M 



(2.11) 



M jk ■ 



In order to find the Euclidean gauge action, we must keep track of the temporal vs spatial 
directions. Consider the anisotropic Minkowski lattice spacings defined by: 



x+0 



M 



x+j 



dz 3 



^ = ar/ar + 0(a) 



where ^ is the bare anisotropjfl of the lattice. The corresponding Euclidean lattice 
spacings are given by: 



,M 



£0 = 



The anisotropic gauge action becomes 



7^Tr(l - +^^ReTr(l - 

.j<k ^ 3 



J2 E ^aJ-TriF^) + 0(a 5 ) 



X fJ,,U 

-Sg 



where the anisotropic Euclidean gauge action is given by |Kla98a] : 



Sr, 



N 



E 



V ^ReTr(l - U jk (x)) + V£°ReTr(l - U jA (x)) 



.j<k 



(2.12) 



1 The bare anisotropy parameter £q in the action is note the same as a s /a T due to quantum corrections. 
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2.4 The Euclidean fermion lattice action 

For clarity, we will show all sums explicitly in this section and begin by describing the 
fermion action on an isotropic lattice. On the lattice, we will replace first and second 
covariant derivatives by symmetric differences using gauge links as parallel transporters: 

<9 M -> -(V^)^ = — {U ll (x)5 x+fl ,y - U ] (x - jj)8 x ^ y ) 

— (A^)^ ee — {U^(x)5 x+fl , y + Uj t (x - p)6 x ..fi tV - 25 X:V ) 

We use symmetric differences in the action in order to guarantee a Hamiltonian which 
satisfies reflection positivity: 

H T ^ T H. 

Minkowski time evolution e tHt O(0)e~ lHt requires H* = H, but Euclidean time evolution 
e Hl ~O(0)e~ HT requires reflection positivity. 

Making the appropriate replacements, we may write the fermionic part of the lattice action 

as 

S F = J d 4 xip(x) ^y^jpdp + rnj ip(x), 

= J> 3/2 ?(*)) ( Y^(y») x , y + (am)8 x A (a 3 ^(y)), 

= ( Y^( v ^ x >y + am5x >y) 



where we have defined the dimensionless quark fields ij) x = a 3 / 2 ip(x) and if) x = a 3 ^ 2 ^p(x) and 
have grouped a factor of the lattice spacing a with the bare mass m to make the 
dimensionless quantity am. 



2.4.1 The Wilson term 



The naive replacement of the derivative by the symmetric difference in the above leads to a 
serious lattice artifact known as fermion doubling [Rot 9 7] in which there are 15 additional 
massless fermion 'doublers' which contribute to the calculations. These doublers can be 
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made harmless by introducing the Wilson term, an 0(a) term (irrelevant in the naive 
continuum limit) into the action which gives the doublers a mass of 1/a: 

S F -> S F - a- ^2 a^i^d^d^ipiy) 

x,y,n 

= 

Q = am + J^(7^V M - -A M ) 
Qx, y = {am + A)S x> y - [i 1 ~ l^ u ^{ x ) s x+fi,y + (i + l^)Ul{x - ^S^y] 



2.4.2 Anisotropic lattices 

We would like to use large spatial volumes in order to mitigate finite volume effects in our 
calculations. This requires a coarse spatial mesh in order to have a manageable number of 
spatial lattice sites. On the other hand, our temporal correlation functions decay 
exponentially with the baryon mass, motivating us to work with a fine mesh in the 
temporal direction. Accordingly, we will work with an anisotropic 

lattice [Kla98a, IMP991 Kla99] with spatial lattice spacing a s , and temporal lattice spacing 
a T . The anisotropic gauge action is given by Eqn. 12.121 Repeating the derivation above for 
an anisotropic fermionic lattice action is straightforward if we introduce the elementary 
lattice vector a M = (a s , a s , a s , a T ). We now absorb factors of a^J 2 into tp and ip, and 
associate a factor of a T with m. In addition, we put in an explicit factor of the bare speed 
of light v which takes into account the fermion anisotropy effects. This gives us 

Sf = 1PQ1P 

Q = a T m+-J2(ljVj " + (74V4 - -A 4 ) 

Qx,y = (a T m + — + l)O x ,y 

4o 



(1 - r )j)Uj(x)5 x+ iy + (1 + 1] )U]{x - j)5 x . 

240 ^— ' 



J-U 



J 



l - [(1 - i,)U A (x)5 x+ly + (1 + lA )U\{x - l)5 x _ ly ] (2.13) 
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2.4.3 Setting the scale 

In order to get back to energy units (i.e. MeV and fm) from lattice units (a T = 1), we will 
relate a measured quantity such as the heavy quarkonium string tension (which is not very 
sensitive to mass tuning) to an experimentally measured quantity. 

The heavy quark string tension for large separation R is given experimentally by: 

V(R) &V + aR, 

where a ~ (465 MeV) 2 . On the lattice, we can measure the static quark- ant iquark 
potential and fit to the dimensionless form: 

a T V(R/a s ) = (a T V ) + (a T a s a)(R/a s ) 
= (a T V ) + ^a 2 T a)(R/a s ) 

= (a T V ) + $a){R/a a ). 

The anisotropy can be determined by comparing the behavior of the static quark-antiquark 
potential as a function of both spatial separation (regular) and temporal separation 
(sideways). 

Once the anisotropy has been determined, the string tension can be used to find the 
temporal (or spatial) lattice spacing in units of MeV -1 . This fixes the energy scale of all 
quantities in the simulation [EHK98a, LM93J. 



2.5 Tuning the lattice action 

The general steps involved in tuning the lattice paramters typically involve: 

• Tune the bare anisotropy £o and compare the regular and sideways static quark 
potential to determine the renormalized anisotropy £ = a s /a T [Kla98aJ. 

• Tune the bare coupling (3 and use the string tension (y/a = 465 MeV) to determine 
the temporal lattice spacing a T or spatial lattice spacing a s |EHK98aJ. 

• Tune the bare fermion anisotropy v (the 'bare velocity of light') and measure the 
pion dispersion relation, which will converge to E 2 = m 2 + p 2 as the fermion 
anisotropy approaches £ |Kla99[ ISLLL06] . 

• Tune the quark masses and measure the spectrum. As quark masses decrease, 
spectral quantities should scale according to chiral perturbation theory. 
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For the current work, we use unrealistically high quark masses in order to have rapid 
convergence of our propagator inversions. The pion mass can be estimated by a 
back-of-the-envelope calculation: 



where we have rounded to the most significant digit due to the crudeness of our estimate. 

2.5.1 Lattice parameters for this study 

Our study uses the anisotropic Wilson action with the following parameters. 

• Lattice size: N s = 12, N T = 48 

• fo = 2.464 <-> f « 3.0 

• (3 = 6.1 <-> a s pa 0.1 fm 

• a T m ud = -0.305 «-> « 700 MeV 

• i/ = 0.902 

The configurations used for this work were quenched (discussed in the next chapter) and 
the spatial volume considered was relatively small (~ 1.2 fm). This was acceptable because 
the purpose of the present work was to develop and test a new operator design 
methodology, specifically with the goal of designing and tuning good nucleon operators to 
be used in later 'production' runs to extract the low-lying baryon spectrum in the nucleon 
sector. It is expected that the operator design methodology described in this work will 
remain relevant as we work to larger spatial volumes, lower quark masses, and finer lattice 
spacings. Production runs will use the anisotropic Wilson action with dynamical u, d, and 
s quarks and a clover improvement term [A DL + 95l IAKL97t IKla98bj IEHR03] . 




(from sideways potential tuning) 

(from string tension) 
(from pion effective mass) 

(he pa 200 Mev-fm) 



a s pa 0.1 fm 
a T m n = 0.1125(26) 

pa 700 Mev 
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2.6 The spectral representation of correlation 
functions 



The Euclidean formulation of a quantum field theory provides a simple way to calculate 
vacuum expectation values [MM94] . Consider 



Tr(e- HT A) 
Ti(e- Hr ) 



J2(k\e- Hr A\k) 



k=0 



1=0 

oo 

Y,(k\A\k)e- E > 



k=0 



1=0 

™ (n\A\n) 

To see how this Hamiltonian approach connects to the functional integral approach with a 
lattice regulator, we will use a complete set of field eigenstates {\<fi)} in the trace, rather 
than the complete set of energy eigenstates For clarity we will consider the case of a 

complex scalar field operator $ for the rest of this section. Let the integration measure 
over field configurations be denoted by dfi(4>). Then for a fixed time slice we have: 

$\<t>) = 



10X01 = / #*#e-^|0)<0| 5 



where the factor of ^ arises due to the definition of the field eigenstates. This will be 
illustrated in the case of fermions in Subsection 12.6.11 A rigorous treatment for both 
bosons and fermions is given in [Bro92j . Consider the following partition function with the 
trace taken with respect to a complete basis of field configurations: 

Z = Tr(e~" T ), 

d/i(0 o ) (0o|e" HT |0 o ), 
(inserting N T — 1 complete sets of states, with T = N T a T ), 

dli{<t>N T -xW{<l>N T -2) ■ • -d/i(0o) (<Po\e~ HaT \<pN T -i)(<pN T -i\e- HaT \0N T -2) ■ ■ ■ (0i|e- Ha -|0o), 
dfi(4> NT ^) ■ • -g^(0o) {(j)N T \e- HaT \(j>N T -i) • ■ ■ (0i|e-^|0 o ). 

i>N T =<f>o 
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It can be shown [Bro92j that: 



-Hot 



\4> T ) = exp{0; +1 (f)0 T+1 (x) - a T L(<f) t+1 (x),(f) t (x))}. 



(2.14) 



The quantity on the left is the matrix element of the transfer operator e~ T between two 
states in the Hilbert space. The quantity on the right is a function of two different field 
configurations <fi T+ i(x) and <fi T {x). We have shown the spatial indices as a reminder that we 
are dealing with the field degrees of freedom defined over all of space (the Hamiltonian is 
time-independent). If we regularize the individual integrals over the (f) T by placing the 
system in a finite spatial volume with periodic boundary conditions and by discretizing the 
spatial coordinates we have: 



The choice of periodic boundary conditions (p.b.c.) for the spatial volume eliminates edge 
effects, but does does not eliminate finite volume effects due to images. The periodic 
boundary conditions in the 'temporal' direction (r) arise naturally as a consequence of the 
trace. It is not necessary to use a trace to extract the spectrum, but the trace has the 
useful property of temporal translation invariance. The anti-periodic temporal fermion 
boundary conditions to be used with the trace method are discussed in Subsection 12.6.11 

If we identify the lattice action S[(f>] = J2 T a rL(4>T+i,(pT), we arrive at the fundamental 
relation between the trace of the operator e~ HT and the (well-defined) functional integral 
over field configurations <fi T (x) = <fi{x) on our Euclidean space-time lattice (identifying 
x A = r): 



For finite lattice spacings a T} a s and a given lattice action S, the lattice Hamiltonian H 
defined via Eqn. 12.151 will differ from the continuum Hamiltonian of the theory by terms of 
order a T and a s . The goal of action improvement programs is to reduce these discretization 
errors by adding offsetting 0(a) terms to the lattice action S. 

We can access the spectrum of the theory by inserting an arbitrary creation operator O at 
'time slice' and an arbitrary annihilation operator O at 'time slice' r into the functional 





p.b.c. 



(2.15) 
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integral expression for Z: 



(C(r)C(O)) = \ ( P0C(r)O(O)e- 5M (2.16) 



Z 



p.b.c. 



Tr (e- H ^Oe~ HT 0) 



Tr (e- HT ) 





e -H(T- 




- Ht o\ 


n) 




e~ HT \ 


m) 



En e- EnT (n\e +HT Oe- HT 0\n) 



(Q\e +HT Oe- HT 0\Q) (2.17) 

oo 

^2(n\Oe- HT \k)(k\o\n) 



k=0 

oo 



J2\(k\0\m 2 e- EkT (2.18) 



fc=i 



where in 12.171 we have taken the limit T — > oo (large temporal extent of our lattice). If we 
require our creation operators O to have a vacuum expectation value of zero 
((0|(5|f2) = 0), then we can extract the first (non-zero) energy state of the theory by fitting 
to the asymptotic form of Eqn. 12.181 

The form of 12.171 looks like a Heisenberg time-evolution equation if we replace t — > —ir 
(giving diffusion dynamics, not wave dynamics). This was to be expected because we 
constructed our Euclidean theory such that the correlation functions would be related to 
the Minkowski Green's functions analytically continued to imaginary time. 



2.6.1 Temporal boundary conditions for fermions 

There is some subtlety involved with taking the trace using eigenstates of fermion fields 
due to their Fermi-Dirac statistics |Bro92j . For clarity in this section, we consider the 
mathematics of a single fermion field carrying no indices. 

The fermion fields in the functional integral are Grassmann (anticommuting) numbers, not 
c-numbers (which commute). Consider the case of two independent Grassmann variables: z 
and z* (not related by complex conjugation). These objects satisfy: 

zz = —z z, z = z =0. 

Consequently, any polynomial in these variables has a finite (exact) power series expansion: 



P(z, z*) = po+ Piz + P2Z* + Vvizz* . 
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Specifically: 

e 2 = 1 + z. 

We can define integration over Grassmann variables formally by 

dz = I dzz=\ [ dz* = I dz*z* = 1 



with the requirement that j dz* dz = — J dz dz* for consistency. Note that the above 
definition defines a translationally invariant integral 

J dz (z + d) = J dz z 

where d is any c-number or Grassmann number. 

The pair of fermionic field operators \l/ and ^ satisfy the fermion anticommutation 
relations: 

V]>2 _ ^f2 _ q _|_ \J/txJ/ _ I 

These relations are sufficient to show that these operators act upon a two-dimensional 
Hilbert space spanned by the states |0) and |1). The behavior of the operators \l/ and ^ 
on this space is completely specified by 



and 



*|0) = 


(0|*t = 




*t]l) = 


= 




(0|0) = 


(111) = 


1, 


(0|1) = 


(1|0) = 





M/t|0) = |l), 


(0|tf = ■ 


(11, 


*|1> = |0>, 


(ll^t = 


(0|. 



In order to take a trace over field configurations, we will need the eigenstates of the field 
operators. These are given by: 

\ z ) = e ^ z \0) (z\ = (0\e z ^ 
where the Grassmann numbers z and z* anticommute with the operators \I/ and ^ 

\ff Z = -zV, Z^ j = -^Z, 

\$J Z * = - Z *^f, z*^ = -&z*. 

Using the rules we have just defined, we can verify that 

ty\z) = z\z), (z\^ j = (z\z*. 
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Notice that (by construction) this looks like complex conjugation. To define the measure 
over the set of field eigenstates, we need 

( Zl \ Z2 ) = (0|e^*e* t22 |0) 

= (0|(1 + ^)(1 + ^ 2 )|0) 

= (0|0) + ^(0|^|0) + (0(^10) + ^2(01^10) 
= l + z *z 2 = e z * Z2 

If we define the measure over coherent states dfj,(z) by 

dfx(z) = dz*dz e~ z z 

then we have the resolution of the identity in terms of a complete set of coherent field 
states: 

d[i(z) \z)(z\ = 1. 



Checking this: 

J dfi{z){ Zl \z){z\z 2 ) = J dz*dze- z * z e z ^ z e z ' Z2 , 

= J dz*dz (1 - z*z){\ + z{z){\ + z*z 2 ), 

= J dz* z* dz z (1 + z*z 2 ), 
= e z > 2 = ( Zl \z 2 ). 

Alternatively: 

J dfi(z)\z)(z\ = J dz*dze- z * z (l + ^z)\0)(0\(l + z*ij) 
= |0)(0| + |1)(1|. 

Let F(ty, \l/t) be an analytic function of an even number of fermion operators SE^ and an 
even number of antifermion operators The trace of F is: 

Tr(F) = J dn(z) (~z\F\z) 

= J dz*dz (1 - 2*z)(0|(l - z**)F(l + &z)\0) 
= (0|F|0) + (1|F|1) 

Notice that the trace must be taken between (—z\ ■ ■ ■ \z), not (z\ ■ • ■ \z). In the functional 
integral, this translates to using periodic temporal (rr 4 ) boundary conditions for the 
bosonic fields (i.e. gauge links) but anti-periodic temporal boundary conditions for the 
fermionic fields (i.e. quark fields). We use periodic spatial boundary conditions for both 
bosonic and fermionic fields. These boundary conditions will be implied when we write 
'p.b.c.' on our functional integrals. 
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2.7 Effective mass plots 

Returning to the spectral representation of the Euclidean correlation function: 

oo 

C(r) = (O(t)O(0)) T =>°° \ck\ 2 e~ Ek \ c k = (k\0\Q), (2.19) 

k=i 

we can get an idea of how well the temporal extent of our lattice is approximating the 
T — > oo limit by examining how calculated physical quantities vary with the number of 
lattice sites in the temporal direction for a fixed temporal lattice spacing a T . 

We wish to extract as many of the energy levels Ek as possible from the exponentially 
decaying signal in Eqn. 12.191 It will be a recurring theme in this work to examine plots of 
effective mass functions to judge the quality of our signals and to tune our parameters. 
The effective mass function is defined by: 

cl t M(t) = In (C(t)/C(t + Or)) ■ (2.20) 

This definition of the effective mass is not unique, but it is one of the simplest, involving 
the value of the correlation function at only two neighboring values of time. 

For large values of r in the T — > oo limit, the effective mass may be expanded as: 

C(r) 



a T M(r) = In 



In 



C(r + a T ) 





Cl 


2 e -E 1T + 


c 2 


2 e -E 2 r + . . . 


|ci 


2g-El (Y+Clr) _|_ 


c 2 


2 e -E 2 (r+a T ) _)_... 



= a T E 1 {l + O {\c 2 \ 2 /\ Cl \ 2 e-^~ E ^^)) 

and thus a T M(r) approaches a plateau of a T E\ for sufficiently large r/a T . We emphasize 
that effective mass plots are simply visualization tools; final fits are best made to the 
correlation matrix elements directly. 

Because the lattice is periodic (anti-periodic) in the temporal direction for bosons 
(fermions), we can expect significant contamination from the backward propagating states 
for large correlation lengths r. In the case of meson operators, the backward propagating 
state contamination comes from a time-reversed mirror image of the meson itself. We can 
explicitly take this contamination into account by defining the effective mass o 1 -A{/ cos i 1 (t) as 
the solution to: 

ln[cosh((a r M cosh )(r/a T ))] = ln[cosh((a r M cosh )(r/a T ))] + const. 



An example of both types of effective mass plot is given in Figure 12.11 which gives the pion 
effective mass as measured on our configurations. 



50 



0.4 



0.3 



0.2 



0.1 



0- 



(3 = 6.1 ^ = 3 12x48 
a T m ud = -0.305 N cfg= 100 



at 



m 



m 



m 



am =0.1125(26) 



[] 



m m m 



m 



m in 



m 



32 



m 



m 



_l I I L_ 



_l I I L_ 



_l I I L_ 



_l I I L_ 







10 



15 



20 



x / a 



Figure 2.1: The pion effective mass plots. The standard effective mass a T M(r) (red) fails to 
plateau due to backward channel contamination. This contamination is explicitly treated by 
the meson effective mass function a T M cos h(T) (blue). A cosh fit to the correlation function 
yields a pion mass in lattice units of a-rin^ = 0.1125(26). 



In the case of baryon operators, the backward propagating state contamination comes from 
states in the opposite-parity channel, which will generally have different masses due to 
spontaneous chiral symmetry breaking. As a result, we are unable to model the correlation 
function using a simple cosh form. In this work, we do not correct for this contamination, 
but limit our effective mass plots to the first or twenty values of r/a T (we are working on a 
lattice with N T = 48). Because our baryon states are much heavier than mesons, the 
backward propagating states decay significantly by the time their signal reaches early 
values of r/a T . This is treated in detail in Chapters and [SJ 
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2.8 Extracting excited states 

In order to extract higher states, we may try to fit multiple exponentials to the decaying 
correlation function. In practice, the correlation function is estimated using numerical 
techniques (c.f. Chapter [3]), which results in an uncertainty associated with each data 
point. The method of fitting many exponentials to a numerically estimated correlation 
function is ineffective due to the rapid decay of the signal to noise ratio. 

We introduce a method for extracting multiple excited states which is based on the 
approach of Liischer and Wolff [LW90J. In this approach we utilize a set of n 'candidate' 
operators {Oi, O2, • • • , O n } and construct linear combinations 

n 

®i = O a v ai , 
0=1 

with associated correlation functions: 

Ci(r) = (6^)6,(0)), 

= vjC(r)vi, 

where we have introduced the correlation matrix C a b(r) = (O tt (r)C?&(0)) and the fixed 
coefficient vector v\. C(r) consists of correlation functions between every possible pair of 
candidate operators in our original set, and therefore contains more information about the 
spectrum than the correlation function of a single operator. Note that all of the Cj(r) and 
C a a( r ) (diagonal elements) are real and positive for all r ^0. Our task is to extract as 
much information as possible about the spectrum from the correlation matrix C(r) by 
finding suitable coefficient vectors {v^}. 

2.8.1 The variational method 

Let the state \vi) be defined by: 

\ Vi ) = e- H ^ 2 e i \n) = e- H ^ 2 J20a\n)v ai 

a 

(Vi\ = (Q\Q l e- H ^ 2 = J2<(^Pae- HT0/2 

a 

for some reference time r > 0. In inclusion of a reference time in the normalization of \vi) 
is necessary to ensure that we never have to evaluate a correlation function at zero time 
separation: 

( Vi \vi) = (fiie.e-^e^) 

= a(T ) 

= vjC(T )Vi 
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(v 


F 


v) 


(v 


v) 



C(0) will generally involve products of multiple fields evaluated at a single space-time 
point. Such expressions lead to subtle difficulties which arise during the regularization 
process [Wei95j. Specifically, C(0) will often receive nonzero contributions from 'Schwinger 
terms' |Sch59j . 

In the space of all possible states {|fj)} parameterized by the vectors V{, we would like to 
find states which are as close to the stationary states {\k)} as possible. The variational 
method [CT77] is a powerful technique that is used to diagonalize some Hermitian operator 
F within some subspace spanned by a set of parameterized trial states {|t>)}. To illustrate 
this method, let f(v) be the expectation value of F in the state \v)\ 

f(v) -- 

(v\F\v) = f(v)(v\v). 

Taking the first-order variation with respect to the parameter v gives 

(Sf(v))(v\v)+f(v)(5v\v)+f(v)(v\5v) 
(5f(v))(v\v) 

where we have used the Hermitian property F* = F. As long as (v\v) 7^ 0, we see that 

Sf(v) = Q <- F\v) = f(v)\v). 

The variational method diagonalizes the operator F restricted to the subspace spanned by 
the If any eigenstates of the unrestricted operator F are within the subspace, this 

method will find them. Otherwise, the method will produce the 'closest states' to the 
eigenstates of the full operator F it can find in the subspace. 

For our purposes, we will vary the parameters Vi until the following expression is at an 
extremum: 



(Sv\F\v) + (v\F\5v) 
(6v\(F-f(v))\v) + (v\(F*-r(v)\6v) 



(Vi 


e -H( 


T— To) 


Vi) 


(Vi 


Vi) 



v\C{r Q )vi 



The infinitesimal change 



gives 



vlC(r)v t 

vjC(To)Vi_ 

(yi + dViYCtfivi + dVi) 
<)r< C(f)Vi- 



{vi + 5vi^C(T)(vi + Svi) v\C(r)vi 



(vi + 5viyC(T )(vi + 5v. 
vjC(r)v 



vlC(T )Vi 



vlC(r )vi 

v\C(T Q )Vi 



Vi + 8v^C{T Q )(vi + 5vi 



C(r )vi 



+ (c.c.) + 0(Svt 



0, (2.21) 
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where (c.c.) stands for the complex conjugate of the first term, and we have used the fact 
that Ct(r) = C(r). 

Equation (I2.2ip is satisfied for all infinitesimal variations 5vi if: 

C{r)v i = X i (T,T )C(T )vi (2.22) 

The eigenvalues Aj(r, To) (ordered by increasing value at each r) of this generalized 
eigenvalue problem are called the principal correlation functions. 

To find the asymptotic behavior of the principal correlation functions expressed in terms of 
the desired energies E k , consider the quantity v\C(j)vj. From the orthogonality property 
of the generalized Hermitian eigenvalue problem, we can always normalize our eigenvectors 
Vi such that: 

vlC(T )Vj = (Vi\Vj) = 5ij, 

which gives: 

v l C ( T ) v 3 = ^j( T > T o)vlC(T )v j} 
(wile-^^l^) = i^j. 

We can write any state in the Hilbert space as a linear combination of the energy basis 
states 

\ V j) = \ k ) a ^V 
k 

with the associated time evolution (or relaxation) given by 

e -H(T-T )\ v .} = \k)a kj e- Ek{T - T0 \ 

k 

Here we see that \vj) evolves to a state which is orthogonal to \vi). Because \vi) has no 
overlap with the large-time relaxation state of \vj), it is forced to relax to a different 
excited state. Thus, the n principal correlation functions A;(t) behave as: 

lim A fc (r, r ) = e - Ek{T ~ T0 \l + 0{e- Ak(T - To) ) (2.23) 
A fc = min \E t - E k \ 

l^k 

where the E k are the energies of the first n states accessible from the vacuum by the action 
of our creation operators O. An energy state \k) is accessible if there exists some 
a G 1, • • ■ , n such that 

(k\o a \n)^o. 

The contamination arises due to the fact we are restricted to the subspace of states 
spanned by {O a |f2}}. 
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Looking at Eqn. 12.231 we may define the principal effective mass functions: 

Afc(r,r ) 



a T M k (r) 
lim a r Mk(r) 



m , 

_X k (T + a T ,T ) 

a T Ek 



We may turn the generalized Hermitian eigenvalue problem (12.221) into a regular Hermitian 
eigenvalue problem: 



C(rM 

c-^i^c^c- 1 / 2 ^)^ 1 / 2 ^ 

[C-^^C^C- 1 ' 2 ^)] Ui (r,r ) 



Xi(T,T )C(T )Vi, 
\i(T,T )C(T )Vi, 



(2.24) 



where we assume a positive-definite correlation matrix at r = r and have added an explicit 
r dependence to the eigenvector solutions Uiij, r ) to remind the reader that fluctuations in 
the Monte Carlo estimates of the correlation matrix elements at different values of r will 
result in corresponding fluctuations in the solution of Eqn. 12.241 

In practice, degeneracies or near-degeneracies in the energy levels combine with numerical 
uncertainties to introduce ambiguities into the estimate of a consistent set of 
fixed-coefficient vectors Vi and the consistent ordering of the principal correlation functions 
Aj across different values of r. In order to avoid these ambiguities and to fix a basis in 
degenerate subspaces, we will adopt the approach of solving for the fixed coefficient vectors 
Vi at a single time r*: 

v t = C-^ 2 (r )u(r*,r ), 

where t* is chosen to be as small as possible (to avoid noise) but large enough to ensure 
that the estimates of the Vi have stabilized. The diagonal elements of the rotated 
correlation matrix C(r) are called fixed- coefficient correlation functions, and are defined by 

Cfcfc(r) = vlC(r)v k , 

and can be used to define fixed- coefficient effective mass functions: 



a T M k (r) 



In 



Ckkij) 



C k k(r + Or) 



In 



v\C{t + a T )v k 



The fixed-coefficient correlation functions Ckk{ T ) differ from the principal correlation 
functions Afc(r, r ) because we are only solving the eigenvalue problem once, at r*. Thus we 
are only guaranteed orthogonality at r*: 



Cij(T*) 







* 7^ 3- 
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The lack of orthogonality for r > r*, no matter how small, means that our fixed-coefficient 
correlation functions will eventually relax to the lowest excited state 

Ckk{r) T -^ |ci| 2 exp(-£ir). 

In Chapter [SJ we will seek (and find) a range of r values over which orthogonality holds to 
good approximation. This will allow us to extract the values of the excited states using 
exponential fits to the fixed-coefficient correlation functions. 

2.8.2 Discussion 

The variational method allows the physics to tell us which combination of the candidate 
operators best interpolates for the stationary states of the Hamiltonian. Because the 
coefficients {v a i} are time-independent, we may use them to infer something about the 
physical structure of the state \k) from the way in which the operators O a combine to form 

e fc . 

Our need to perform the inversion and square root in Eqn. 12.241 underscores the importance 
of choosing candidate operators which give a correlation matrix with a good condition 
number. We will keep this in mind in Chapter [7] when we prune our large operator set 
down to manageable size. 

Our task is to construct a 'good' set of candidate operators O a which we believe capture a 
representative set of features of the baryon spectrum we seek (in order to construct a 
'good' trial subspace for the variational method). 



Chapter 3 

The Monte Carlo method 



Our goal is to compute the low-lying spectrum of lattice QCD. We do this by calculating 
the elements of Euclidean space-time correlation matrices, each of which is a ratio of 
functional integrals: 



C ab {r) 



where 



Q*,v[u] 



s G [u] 



(O a (r)O b (0)), 

\ [ VUV^V^O a (r)O b {Q)e~ SF ^^- s ^ u \ 

Z ip.b.c. 



VU Vtfj Vtfj e -SF®,4>,v\-So\u] ? 

p.b.c. 

4Q[U]il>, 

(a T m + — + l)d x> y 



1 - l4)U A {x)5 x+ly + (1 + ^M(x - A)5 X _^ 



I 
N 



E 



V iReTr(l - U jk {x)) + V^°ReTr(l - U j4 (x)) 



(3.1) 
(3.2) 



(3.3) 

(3.4) 
(3.5) 



(3.6) 
(3.7) 



Before we turn to the design of the operators O a , we briefly discuss how we will evaluate 
the ratio of functional integrals defined by Eqns. 13.21 and 13.31 
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3.1 The need for the Monte Carlo method 

These integrals are regularized by the space-time lattice, but still consist of an 
unmanageable number of integration variables. The Feynman diagram method is not 
applicable to QCD functional integrals because expansions parameterized by the coupling 
constant do not converge at the energy scales in which we are interested. 

We turn instead to the Monte Carlo method, a powerful technique which requires us to 
formulate the problem on a computer. Our quark fields satisfy Fermi-Dirac statistics, and 
are therefore represented by anticommuting Grassmann variables which are not suitable for 
direct numerical integration on a computer. Fortunately, the quark fields appear 
quadratically in the action and in the operator pairs O(t)O(0), which means that the 
Gaussian integral over the fermion fields is symbolically integrable. Therefore, we will 
integrate out the Grassmann fields by hand before we translate the functional integral to 
computer language. 



3.2 Integrating the quark fields 



We wish to perform integrals of the form 

V^jVijij Xl ij X2 ■ ■■i) Xn _j\) Xn exp(-^2^ x Q X! y[U]ip y ) (3.8) 
5 5 5 5 



ZfI p . 



bc x,y 



5J X2 5J Xl $Jx n 5J Xn _ 1 



\n(Z F [J, J])\ J>J=0 



where we have introduced the fermion generating functional Z F [J, J}: 

Z F [J,J} = [ V^Vi;e^l-J2^xQxAU]^y + Yl^ + ^ J ^ ) > ( 3 ' 9 ) 

J P- hc - \ x,y x J 

which will allows us to insert arbitrary factors of ip x ip y into the functional integral by 
taking functional derivatives of Z F [J, J] with respect to the source fields J y and J x and 
then setting J y = 0, J x = 0, which eliminates the source terms and gives Z F = Z F [0 : 0]. 

The sources J X) J x and the functional derivatives ^j-, jj- obey the Grassmann algebra just 
like the quark variables ip x and ip x , and satisfy: 

5 5 5 5 

jjjv = s *,v jj x J y = Jjx Jy = 6x > y 5lx Jy = 

We can integrate Z F [J, J] by shifting the quark field variables to complete the square in 
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the exponent: 



y a = c,.-J !l Q, i -.r 

"fix = ^x-QZl[U]Jy 



z F [j,J] = [ vyvtf expl-J2^ x Q*Mi>y + J2 J *Q*iW J y 

■'■» hA - V x,y x,y 

= / IJ^' d ^ x ' eXP ( ~^2^xQx,y[U]^ V I exp I -JxQx!y[U}Jy I (3.10) 

■y 1 \ t n / \ nr. it / 



x,y / \ x,y 



where we have dropped the primes, written the integration measure explicitly, and have 
used the fact that the product of two Grassmann numbers behaves as a c-number. The 
product of integrations extends over all lattice sites x', and the fermionic boundary 
conditions are applied to expressions such as 5 XjX+ fi within the quark matrix Q[U] and 
inverse quark matrix Q~ 1] \U]. 

The Gaussian integral over Grassmann variables in Eqn. 13.101 can be evaluated [MM94] by 
expanding the integral and keeping only terms which do not vanish when integrated: 

J J] exp [-Y^xQxA^y] = det(Q[E/]) / J] (# x ,Vv) (dWW) = det(Q[U}) 

x' \ x,y / x' 

Thus 

Z F [J, J] = det(Q[V)) exp (j2 J x Q^ y [U]J y 

V x,y 



At this point it is essential to address the correct treatment of the suppressed color, flavor, 
and Dirac indices. The matrix Q[U] carries those indices, as do the variables ip, ip, J and 
J. The quark matrix Q[U] can be thought of as a tensor product of matrices acting in 
position space, color space, flavor space, and spin space: 

• 8x,x+fi acts in position space 

• (U^) a b(x) acts in color space 

• Q[U] is diagonal in flavor space (QCD conserves flavor) 

• {in) a/3 acts in spin space 



Consequently, when we take a matrix determinant or inversion, we must do so with respect 
to all indices. 
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We can use the fermion generating functional to define the full generating functional 

z[j,jy. 

Z[J,J} = I VU det(Q[U])e- SG ^ +jQ ^ J (3.11) 

J p.b.c. 

where now the functional integral is only over gauge links, with the fermionic boundary 
conditions applied to Q[U], Q _1 [i7], and the bosonic boundary conditions applied to Sg[U]. 

Any correlation matrix element of the form in Eqn. 13.21 can be evaluated with the use of 
the generating functional in Eqn. 13. 1 1 L For every fermion pair ipxipy, we apply jj^jj- to 
ln(Z[J, J]). Any gauge link factors are left as they are. We then set the Grassmann source 
fields J and J to zero, and evaluate the remaining regularized and well-defined integral 
over the gauge links U. 

For example, the single quark propagator is: 

(My) 

z 

In the case of two quark/ ant iquark fields, an additional term is generated by the 
derivatives: 

(^ X2 ^J = \\ VU (Q^ X2 Q~l X4 - Q,.!.rX),l,.;) det(Q)e"^ 

The general way to evaluate correlation functions of arbitrary numbers of quark-antiquark 
field pairs is to add up all possible permutations the fields which are in alternating ip-ip 
form (this can be done graphically using Wick's theorem [Bro92j). For each term we can 
read off the inverse quark matrix indices from the indices of the fields, and include a minus 
sign if we performed an odd number of anticommutations to get the quark fields into their 
final ordering. 

In Chapter [5], we will need correlation matrix elements of the form 
The six relevant permutations and signs are given by: 



Permutation 


Sign 


14 25 36 




14 26 35 


+ 


15 24 36 


+ 


15 26 34 




16 24 35 




16 25 34 


+ 



1 

z 



VUQ-i[U)det(Q[U))e- s °W 

p.b.c. 

VUdet{Q[U])e- SG[u] . 



p.b.c. 
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Thus 



p.b.c. 



X'G 



Qx lt X4,Qx2,Xe,Qx3,X5 
Qxi, X$Qx2, X^Qx3, Xq 

- QxUQ 

X2,X6^Z X'j,X4, 
Qx 1 ,X§Qx2,X4,Qx2,,X$ 

+ Qxtx,Qxix,Qxlx A )^{Q)e- SG (3.12) 



The integrand now contains only pure functions of the gauge link variables: Q[Z7], Q 1 [U], 
and S G [U). 



3.3 The Monte Carlo method of integration 



Now that we have a well defined integral and have integrated out the quark degrees of 
freedom, we would like to evaluate expressions such as: 



(f[U\) = I / VU f[U}det(Q[U})e 

\ e -Sa[V] 

(3.13) 



'p.b.c. 

fn Xl Mx)f[U]det(Q[U])e- s °W 



lUx 4l dU,(x)det(Q[U})e- s ^ 



where each dU^x) is the appropriate Haar measure over the SU(3) group [Cre85] . There 
are (2 x x N T x 4) SU (3) integration variables, where N s and iV T are the number of 
lattice sites in the spatial and temporal directions, respectively. SU (3) has eight 
generators, so a brute force integration over our small 12 3 x 48 lattice would consist of a 
2 x 12 3 x 48 x 4 x 8 = 5, 308, 416 dimensional real integral over a compact domain. 

The Monte Carlo method was specifically developed to estimate integrals over a large 
number of degrees of freedom [PFTV92] : 



1 

V 



dfj,(x) f(x) 



v 



where 

V = / dji{x). 



v 



The Monte Carlo method uses the normalized integration measure d/j,(x)/V at each point 
in the integration domain to define the measure pp(x)dx over a probability space F. An 
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ensemble of N configurations {xi,x 2 , ■ ■ ■ ,x n } is generated according to the probability 
measure over F, and the integral is estimated by: 



If 1 - 

- / dll(x) f(x) « — 22 / W = f 
V n=l 



The expectation value of the ensemble average / over F is the same as the expectation 
value of the function evaluated on any single member of the ensemble Xj, and equals the 
desired integral: 

1 - 

= jyE £ 'W = W)] (3-14) 

n=l 

p(x)dfi(x) f(x) (3.15) 
dfx(x)f(x) (3.16) 

v 



V 

1 

V 



From the central limit theorem |Bil86] we know that if iV is sufficiently largqj, then / will 
be normally distributed with: 

Var F [f] oc 1 (3.17) 

The variance of the estimate Varp[/] is a measure of the uncertainty, or expected error, 
caused by using a finite number N of configurations in the Monte Carlo estimate of the 
integral. We will return to the task of estimating Varp[/] in Section I3"31 



3.3.1 Importance sampling 

For our present purposes, we would like to perform high-dimensional integrals of the form 
in Eqn. 13.131 Because the Euclidean gauge action Sq[U] is real and bounded from below, 
we can absorb the exponential damping term exp(— Sg[U]) into our definition of the 
measure. Weighting the probability of generating a configuration with the Boltzmann-type 
factor exp(— Sg[U]) greatly reduces the variance of our estimates by suppressing 
configurations which give exponentially small contributions to the functional integral. The 
fermion determinant det(Q[?7]) is not real for all configurations because we treat i/j and ip 
as independent variables in the functional integral. This unfortunate situation, known as 
the fermion sign problem, makes it difficult to include the fermion determinant into our 
definition of the measure. 

After absorbing the exponential factor into the measure, we can estimate the functional 
integral in Eqn. 13.131 by first generating a ensemble of lattice- wide gauge link configurations 



Finite ensemble error analysis involves subtle issues and is treated in detail in section ET5l 



62 



{{Un(x)}i, ■ ■ ■ , {Un(x)} N } = {t^ 1 - 1 , • • • , C/^} according to the joint probability density: 

e -s a [uW] 



p 



([/(")) = 



/n^^Wdet(g[C/])e-^^] 
and then approximating: 



OT . i^gD ,3,8) 



JV _ 

' ELidet(g[C/M]) 



3.4 Markov updating 

We now turn to the task of generating the gauge configurations U^ n \ Given some 
configuration, we may generate a new configuration by cycling through the lattice and 
proposing random (valid) changes to the individual link variables. The new proposed 
change is accepted or rejected based on some probabilistic rule which considers the 
Boltzmann factor of the configuration that would result if the link change was accepted. 
The process is repeated for several sweeps through the entire lattice in order to produce an 
updated configuration. This is an example of a Markov chain, a class of processes in which 
the current state depends probabilistically on the previous state. We require our Markov 
chains to be ergodic: any state in the configuration space can be reached from any other 
state in a finite number of updating steps. This is necessary to ensure that we are 
adequately exploring the integration domain during our Monte Carlo integration. 

A Markov chain is initialized to some initial state U^°\ It is then updated until it 
'thermalizes.' At this point, all of the statistical properties of the states being generated in 
the chain have stabilized. For example, one can plot a moving average of the plaquette 
value to determine when the Markov chain has thermalized. The plot will drift rapidly at 
first, and then fluctuate around some average value. If the configuration at which this 
begins to happen is denoted U^ m \ then we take the configurations from f/( m+1 ) to lJ( m + N ) 
for our ensemble. U<® -> [/« -> -> ► -> U {m+l) -> ► u {m+N) 

S v ' 

keep 

Examining the moving average of the plaquette is just a heuristic; in practice 
thermalization is a subtle issue requiring careful treatment. Whenever possible, it is 
preferred to initialize a Markov chain with a thermalized configuration from a previous 
instance of the Markov chain. 



3.4.1 Quenching 



Because the fermion determinant det(<5[[/]) is a complex quantity, it is an extremely 
difficult quantity to estimate numerically due to cancellations. There exist tricks, such as 
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the Weingarten pseudo-fermion method Wci85j, which enable us to include the fermion 
determinant in the updating algorithm as part of the integration measure. Such methods 
are computationally expensive because they involve the inversion of a matrix having 
space-time, color, and spin indices. The inversion of such matrices is currently the single 
largest bottleneck in Lattice QCD calculations. 

For this work, we will work in the quenched 'approximation' where we simply set the 
fermion determinant in Eqn. 13.181 to one: 

det(Q[E/]) -> 1. 

This is an unparameterized 'approximation' which breaks the unitarity of the theory. 
Nonetheless, it is expected to to be mostly harmless for tuning studies with unphysically 
heavy quark masses (such as this work) because quenching is expected to affect the levels 
of the spectrum, but not the coupling of the operators to the spectral states. The fermion 
determinant is a highly non-local object containing the complete contribution to the action 
of the fermionic degrees of freedom. From the (non-convergent) Feynman diagram point of 
view, quenching is equivalent to removing all diagrams containing sea-quark loop 
contributions. 



3.4.2 Configurations used in this work 

Although updating is an important and subtle topic, it is not the focus of this work. We 
briefly describe the method used to generate our gauge configurations. 

The fact that Wilson's formulation of lattice QCD uses compact gauge links removes the 
need for gauge fixing (see the discussion in section [T.7.81) . Thus, it is possible to update the 
gauge links freely according only to the probability distribution exp(—Sa[U]). There is no 
need to require that the new configuration adhere to any fixed gauge prescription. 

The gauge configurations used in this work were generated using the Chroma QCD library 
developed at JLab, and use the parameters reported in Subsection 12.5. 1[ That code used 
the Cabibbo-Marinari heatbath algorithm [CM82, KP85] in conjunction with Creutz's 
overrelaxation technique |Cre87j to update the gauge configurations. 

Because the configurations are quenched, the quark masses did not play a role in the 
updating. However, they are used in the propagator calculations, which involve factors of 
inverse quark matrix elements such as Q~y[U]. The operators presented in this work act on 
the nucleon sector, thus only the bare light quark mass a T m u = a T rrid = —0.305 was 
reported in Subsection 12.5.11 

Quenching will not be used in the production runs, and it will be assumed here that our 
operator choices and tuning parameters will not be significantly affected by quenching. In 
this study, we find that our quenched configurations still yield a useful spectrum. Rather, 
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they serve as an indication of the quality of the operators to be used in this ongoing 
research effort. 



3.5 Error analysis 

Monte Carlo integration is a statistical procedure. It only provides stochastic estimates of 
quantities subject to statistical uncertainty. The use of a Markov chain to generate the 
configurations in our ensemble introduces autocorrelations (discussed below). Additionally, 
we will be transforming our data in a complicated manner (e.g. C~ 1 I 2 (tq)) and need to 
quantify the statistical uncertainty for such quantities. 

3.5.1 Estimate variance in the presence of autocorrelation 

We have an ensemble of N (thermalized) configurations {xi} which looks as though it was 
drawn from some probability distribution F but was actually produced in sequential order 
by a Markov chain. The Markov process introduces serial autocorrelations [Ken06J into the 
ensemble which increases the true variance of the estimate Varp[/]. This implies that 
estimates^] Var [/] of the variance of / will underestimate the true variance Var p [/] of /. In 
other words, our error bars will be too small if our data is autocorrelated. 

To reduce clutter in the following, we define 

g(x) = f(x)-Ep[f(x)}, 

N 

9 = ^Y,9M = f ~ E F [fi 

n=X 

which has the properties 

Ep[g{x)) = 0, 
E F [g(x) 2 ) = Var[/(x)], 
E F [g 2 } = Var[/]. 

2 In this chapter, we will denote sample estimates of statistics with a hat. For example, if we use the 
sample mean of some quantity x to estimate x, then we may write x = x. 
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The variance of the ensemble average is 



Vai F [f] = E F [f 

- N N 

^2 Y Y E F [g{x m )g{x n )l 



N 2 

n=l m=l 



I jv 2 N -i N 

T7i ^2 E F[gM 2 } + ^ Y Y E F[g(x m )g(x n )}, 



N 2 Z-,"^"^ TV 2 

n=l n=l m=n+l 



1 / ? N-l N 

\ 71=1 m=n+l 



£fc.fo(*) 2 ] 



(3.19) 



We can define the autocorrelation function Ck[f(x)] for a function /(x) evaluated on an 
ordered series of observations {x±, x 2 , • ■ ■ , x^} by 

Ck[f{x)] E F [{f(x n+k ) - E F [f{x n+k )]){f{x n ) - E F [f(x n )})\ 



y/E F [(f(x n+k ) - E F \f(x n+k )]) 2 ]E F [(f(x n ) - E[f{x n )]f 
E F [(f(x n+k ) - E F [f(x n+k )]){f(x n ) - E F [f{x n )})} 



E F [(f(x)-E F [f(x)]f 



(3.20) 



There is no n index on the left hand side of Eqn. 13.201 because we assume that the process 
generating the configurations has reached a stationary state (i.e. the statistical properties 
of the sequence are invariant under translations in n). 

We may rewrite equation (13.191) in terms of C k [f(x)}: 

1 / 2 JV_1 - 

\ n=l m=n+l 



= j=V*x F \f(x))(l + l^(N-k)C h \f(x)) 

\ k=l j 

1 / 7V_1 2 Af_1 

= -Var F [/(x)] l + 2^C fe --^fcC fc [/(x)] 

V fc=i fc=i / 

where we remind the reader that Varp[/] is the true variance of /. In the case of a Markov 
chain, C k [f(x)) falls exponentially with k, and for sufficiently large iV we have [BMK73J 

Var F [/] -> 1 Var F [/(x)] (1 + 2A[f(x)\) . (3.21) 

where 

oo 

A[f(x)] = J2c k [f(x)] 

k=l 
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is the 'integrated' autocorrelation of the function /. Thus we see that in the presence of 
autocorrelation, the simple estimate of variance: 

understates the true variance Varp[/]. 

A simple unbiased estimate of Var^[/(x)] on an ensemble {x±,X2, • • , xn} is 

i N 

Var[/(x)] = ^(/(x n ) - ff (3.22) 

n=l 

with a similar expression for C}~. In practice, the Markov chain updating method is tuned 
to reduce the autocorrelation in the data to negligible levels. 

3.5.2 The Jackknife estimate of variance 

The variance estimate in Eqn. 13.221 is applicable to functions / evaluated separately on 
each ensemble member. As the ensemble member fluctuates, so does f(x). Eqn. 13.221 is not 
generally applicable to more complicated functions / which involve multiple members of 
the ensemble simultaneously. 

In the following, each x n denotes a collection of simple quantities (such as correlation 
matrix elements) evaluated on the n th gauge link configuration generated by the Markov 
chain. We are interested in a special class of functions / which depend on the all of the 
ensemble members Xi simultaneously through the ensemble average x: 

1 N 

n=l 

In this work, the x correspond to Monte Carlo estimates of correlation functions, and the / 
correspond to the analysis methods applied to these estimates. Consequently, / will 
typically be a complicated function of the ensemble average x, e.g. taking the inverse 
square root of a matrix or solving a Hermitian eigenvalue problem. We would like to know 
how f(x) varies across ensembles, not across the members of the ensembles. 

The typical variance estimate in Eqn. 13.221 is not generally applicable to functions of the 
mean of an ensemble: f(x). Using estimates f(xi) on individual ensemble members in an 
attempt to estimate the variance of f(x) is often a poor approach because the function 
may only be well defined for certain values of its argument. The individual members of the 
ensemble Xi fluctuate far more than does the ensemble average x, and these larger 
fluctuations may result in an undefined f(x n ). For example, if f(x) = ln(x), then f(x n ) 
will be undefined any time x n < 0, even if x > 0. 
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The jackknife estimate of the variance jEfr82] is a non-parametric alternative that uses 
re-sampling to numerically estimate the expected squared fluctuation of / across 
ensembles. The jackknife method is applicable to virtually any function /, but we will 
specialize it here to functions of ensemble averages f(x). 

Let X( n ) be the ensemble average excluding the n th ensemble member x n : 

1 ^> Nx — x n _ ( x n — x \ 

m=l v 7 



Accordingly, we can define the n th jackknife replicate /( n ) by: 

f(n) = /(%))• 

In contrast to f(x n ) which relies on the information in only the n th ensemble member, the 
jackknife replicate /(„) is supported by the information in all ensemble members except the 
n th . The power of this technique is that we can usually choose N large enough to ensure 
that both f(x) and f{x( n )) are well defined. 

The N jackknife replicates /( n ) are highly correlated, but it is reasonable to expect that 
their mean squared fluctuations are proportional to the mean squared fluctuations of f(x). 
The jackknife estimate of the variance of / across ensembles is: 



/ AT _ ]\ N 

Var jack [/(o:)] = {^ r )J2^ ) -IW 2 

^ ' n=l 



1 - 

A') = AT £ An)- ( 3 - 23 ) 



N , 
n=l 



The jackknife error is defined as the square root of the jackknife estimate of the variance. 

We can verify the constant of proportionality (N — 1)/N in Eqn. 13.231 by considering the 
known case f{x) = x: 



N 

Var jack [^] = ( ^— - ) ^2(x {n) - x ( .)) 2 , 



n=l 



N— 1\ A / fx n — x x x 



E 



N ^\ \N 

n=l 



— X 



1 



N(N 



1 r . .. 



v N 
' ' n=\ 
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which is the expected expression for the (uncorrelated) estimate of the variance of x. We 
also emphasize that we use f(x), not /(.), to estimate f(x): 

f=m. 



A generalization of the jackknife is the bootstrap, in which a large number Nb (typically 
0(1024)) of bootstrap ensembles are each created by selecting N members (with 
replacement) from the original ensemble. / is then evaluated on the bootstrap ensembles to 
yield N B bootstrap replicates f^. The bootstrap estimate of the variance of / across 
ensembles is given by: 



. . / 1 \ N * 

Var boot [/(*)] = (j^y) ~ /[•]) 

j N B 



We use the simple jackknife method to estimate the errors on our effective mass plots and 
the more thorough bootstrap method to estimate the uncertainties in our final fit 
parameter estimates. 



Chapter 4 

Baryon operator construction 



4.1 Operator design goals 

The spectrum is calculated by evaluating correlation matrix elements of quantum operators: 

oc 

C ab (r) = (O a (r)O b (0)) = \ck\ 2 e~ EkT , c k = (k\0\Q), (4.1) 

fc=i 

assuming a large temporal extent T of the lattice and that the lattice action satisfies 
reflection positivity |OS73al [OS 75J. 

This chapter describes the design of operators which excite the low-lying baryon spectral 
states. We seek to find operators such that in Eqn. 14. II is as large as possible for the 
low-lying states of interest while being as small as possible for higher-lying contaminating 
states. Because we are using the Monte Carlo method to estimate Eqn. 14.11 we also need to 
improve our operators to reduce the variance in the estimates of C^ij). With these issues 
in mind, our goal is to design a set of operators which 

• are gauge-invariant, 

• have the correct flavor content, 

• facilitate continuum spin identification, 

• rapidly plateau in the effective mass, 

• have minimal noise in the effective mass, 

• are linearly independent, and 

• compose a set of manageable size. 
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This chapter is based on the operator design approach developed by the LHPC and 
described in [BEF + 05a[ |BEF + Q5bl ILMB+06"] . We begin by introducing the basic building 
blocks of our operators: smeared and covariantly-displaced quark fields, and smeared gauge 
links which are combined into gauge- invariant three-quark operators. Next, we discuss the 
group theoretic approach adopted by the LHPC. The gauge-invariant three-quark 
operators are then endowed with the appropriate flavor structure to make our three-quark 
elemental operators. These elemental operators are then combined into full baryon 
operators which transform irreducibly under spatial lattice symmetry operations (rotations 
and reflections). The chapter concludes by showing how the states excited by these 
operators are connected to the continuum J p baryon states. 



4.2 Building blocks 



—A 



On the lattice, we have quark sources ^„ a (x), quark sinks 4>aa( x )i an d gauge links U^ a b(x) 
where A is a flavor index, a, b are color indices, and a is a Dirac spin index. We have 
returned to writing the space-time index x in parenthesis in order to avoid clutter. Note 
that we are still working with the unitless fields defined previously. 



4.2.1 Spatial inversion 



Under spatial inversion /, the field operators transform as: 

U^Jx, t)U\ = 1pp{-Xi T)I/3a Urtaix, t)U\ = lj)p(-X, r)Ip a 

where the 4x4 parity matrices I and / keep track of how the spin components of the fields 
change under spatial inversion. Summation over repeated indices is implied. To find the 
form of I and /, we consider how the fermion terms in the continuum Euclidean action 
transform under the action of spatial inversion: 

d 4 Xtl) a (x,T)tlj a (x,T) -> 



d 4 x ip a (x, rfjjapdjtppix, t) 



d 4 X1p a (x, T^apdilppix, T) 



d 4 xipp(-x, r)ip y (-x, T)lp a I ia 
d 4 x?p (3 (x, r)^ 7 (x, T)Ip a I ya , 
d A x ipx(—&, r)jj al3 djip a (-x, r)I Xa I a( 3 

- I d 4 Xtp x (x,T)Ix a 'Jja!3la!3dj(x,T)4j CT 

d 4 xip x (-x, r)7 4a/3 9 T V' CT (-x, T)hJ„p 

d 4 X lj) x (x, ^IxaJiapIapdi-lp^X, t) 
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To maintain the covariant derivative, we see that: 

Aj(x) - -M-x), U^x) - u;(-x), 
A 4 (x) -> +A 4 (-x), -v C/ M (-x). 

In order for the action to be invariant under spatial inversion, we must have 

Il T = l, I l3 I T = - lv J 7 4/ T = 74 

Using the definitions of the Euclidean Dirac gamma matrices given in Chapter [21 we see 
that I — I = 74, and that the fields transform as: 

ip(x, t) -> V(-^> r )74, VK^> r ) ^(-^> r )74- 



4.2.2 Gauge invar iance 

We require our hadron operators to be gauge- invariant; all operators must transform as 
scalars under local SU(3) color transformations. We may combine three quark sources (or 
sinks) into a color singlet by contracting the free color indices with the totally 
antisymmetric Levi-Civita tensor: 

+ 1 abc = 123,231,312 
e a bc = { -1 abc = 132,213,321 
otherwise. 



As an example, consider the behavior of the local (single-site) operator under a local gauge 
transformation Q(x) G 577(3): 

e ahc il) a (x)^ b (x)i) c (x) -> e abc i) a \x)il) w (x)i) c ,(x)tt ] a , a (x)£l\, b (x)Q} c , c (x) 

= det (fi f {x))ip a , (x)if) v (x)ip c , {x)e a , Vd 
= i/; a ,(x)ip bl (x)ip c ,(x)e a ' bJc ' 



We will always contract our operators with the Levi-Civita tensor, with summation over 
repeated color indices implied. 



4.2.3 Gaussian quark field smearing 



To successfully extract the spectrum, operators which couple strongly to the low-lying 
states of interest and weakly to the high-lying states must be used. Effective 
phenomenological models tell us that we can improve our operators by associating 
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wavefunctions with the valence quarks, rather than working with point quark sources. This 
translates to the use of spatial quark field smearing in our functional integrals. 

Damping out couplings to the short-wavelength, high-momentum modes is the crucial 
feature which any effective smearing prescription [GLM + 89l lASB + 93] must have. Gaussian 
suppression of the high-momentum modes is perhaps the simplest method one can use. 
Since a Gaussian in momentum space remains a Gaussian in coordinate space, we decided 
to employ a gauge- covariant smearing scheme [AKL96] in which the smeared quark field is 
defined at a given site as a Gaussian-weighted average of the surrounding sites on the same 
time slice (see Figure I4TTT) . 





o o o o 



Figure 4.1: A schematic view of Gaussian quark field smearing. The quark operators are 
replaced by 'fatter' versions which better mimic the 'fuzzy' nature of the quark wavefunctions. 



The smeared quark operators are 'fatter', consisting of weighted contributions from nearby 
values with a (real) Gaussian radius of a. We expect that baryon operators using smeared 
quark fields will better represent the spatial structure of the low-lying QCD baryon states, 
and will therefore interpolate better for those resonances. 

On a given time slice, the smeared continuum quark fields are defined as the integral 
transform of the unsmeared fields using the real, rotationally invariant, Gaussian kernel 
Q(r 2 ): 

^(x) = / d\ Q(r 2 )tfj(x + f), Q(r 2 ) = J_ ■ e^ 2 /^ 2 ), 

J (V27TCT 2 ) 3 

where r 2 = (r • r). The first three moments of the real Gaussian kernel are 

J d 3 r Q(r 2 ) = 1, J d 3 r Q(r 2 )r = 0, J d 3 r Q(r 2 )r 2 = a 2 . 

We can transform the integral equation into a multiplicative operator equation by applying 
the Euclidean spatial Fourier transform relations: 
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to the expression for the smeared quark field: 

4>(x) = - J-^ 3 [ d 3 re- r2/i2a2) ^(x + r} 



(V2 



7TO" 



(V2^f J J (2vr 



1 f d 3 r [ ^ e -^^ 2+2 ^/\-^m 



{V2^) 3 J J (2tt 
(Letting z = r '/a + iak) 

' f d 3 z e~ z2 / 2 ) [ e-° 2k2 l 2 e~ Ck -*f{k), k 2 = k ■ k 



2vr) 3 J J J (2tt) 

^(f) = e a2v2/2 i){x) (4.2) 
Using the relation 

= lim (1 + -) 

Eqn. 14.21 may be approximated on the lattice by introducing the quark smearing matrix S: 

f a 2 \ n " 

$ x = J2S x ^ y} S(a,n a ) = (l + — A J (4.3) 

y a ' 

where 

3 

**,v = E( l7k ( i )U» +I/ 2( i -*)U- 2i '.«) 
fe=i 

£ (Uk(x)6 x+ky -6 x , y ), U- k (x) = Ut(x-k) 

fc=±l,±2,±3 

is the three-dimensional gauge-covariant Laplace operatoiQ. It is simple to check that the 
analogous expression for ip in the continuum 

5(f) = e a2v2/2 ^(f) 

is approximated on the lattice by: 

^ = ^S n ^{a). 

There are two quark field smearing parameters in Eqn. 14.31 that we will tune in Chapter [6j 
the real (dimensionless) Gaussian radius a and the integer number of iterations in the 
approximation of the exponential n a . 

The smeared quark fields at each site are linear combinations of the original lattice quark 
fields; we note that the smeared quark fields obey the same transformation laws and 
Grassmann algebra as the unsmeared fields. 

1 We have defined A and a to be dimensionless. The smearing radius in MeV -1 or fm is given by aa s . 
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4.2.4 Gauge link smearing 

We can improve our operators further by replacing the gauge links in our quark smearing 
matrix S with smeared gauge links U^x). 

Let Ck(x) denote the sum of perpendicular spatial staples which begin at lattice site x and 
terminate at the neighboring site x + fi: 

C k (x) EE J2( U ^ X W X +j) U j( X + ^ 

+ U]{x - j)U k {x - 3)U 3 (x - 3 + fc)) . 

We incorporate the spatial information contained in the staple sums {C^} with the original 
links C/[°l ee {Uk} in some manner to produce a smeared gauge link configuration f/M (the 
superscript refers to the smearing iteration, and has nothing to do with the Markov chain). 
This process is then iterated n p times to produce the final spatially smeared link 
configuration U: 

jjio] _> jjK ^ u {2] -» ► U [np] EE U 

The APE smearing jA + 87j iteration rule is commonly used: 

uf 1] =V8U(s) (yF(x)+pCF(xj) 

where p is some real weight, and Vsu(z) is a projection operator into the gauge group. The 
projection is needed because the sum of two SU(3) matrices is not necessarily an SU(3) 
matrix. 

To avoid the abrupt projection back into the gauge group needed by this prescription, we 
decided instead to use the analytic and computationally efficient gauge link smearing 
scheme known as stout-lin^ smearing |MP04] defined by 

U[ n+l \x) = e W (ipQW(x))U [ k n] (x), (4.4) 

= \ - - ^ (nl(«) - ( 4 -5) 

^k( x ) = Ck( x )Ul(x) (no summation over k) (4-6) 

where N = 3 for SU(3). The two parameters to tune in this smearing procedure are the 
real staple weight p and the integer number of iterations n p . 

To make use of stout links, we replace the quark smearing matrix S in Eqn. 14.31 with the 
complete smearing matrix S: 

( a 2 ~\ na 

ipx = ^2 Sf.v'&v' S ( a ' n a,P,n p ) = ( 1 + - — A J (4.7) 

y a 

2 The term 'stout-link' was coined by Morningstar and Peardon in a Dublin pub as a nod to Guinness 
beer. Officially, 'stout' refers to "their thick-bodied nature from the large brew of paths used in their 
formulation" |MP04j 
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where 

fe=±l,±2,±3 

The quark field and gauge link smearing schemes preserve the gauge invariance of hadron 
operators. In Chapter El we systematically tune the smearing parameters a,n a ,p, and n p 
to give an optimal effective mass signal for our operators. 

4.2.5 Radial structure through covariant displacement 

The need for extended three-quark operators to capture both the radial and orbital 
structure of baryons has been emphasized and described in Ref. [BEF + 05a] and 
Ref. |BEF+05b| . We may add radial structure to our operators without breaking gauge 
invariance by using the covariant displacement operators D^\x) which consist of p 
smeared parallel transporters in the j th spatial direction which 'carry the color' back to the 
reference site x from the displaced site x + pj: 

bf\x) = U 3 {x)U 3 {x + j) ■ ■ ■ U 3 (x + (p - l)j), j = ±1, ±2, ±3. 

For economy of notation in the upcoming expressions, we also define the additional 
zero-displacement operator: 

D<?\x) = l 

A quark sink displaced from x to x + pj is represented by 

Df{x)^{x+pJ) = [Dfi!]{x), 
and a quark source displaced from x to x + pj is represented by 

^(x + pj)dW(x) = $i)f\x) 

For simplicity, the displaced quarks in a baryon operator will always be displaced by an 
equal amount p. In general, operators corresponding to different values of p may be 
combined (via the variational method) to capture more complex radial structure. In this 
work, we restrict our consideration to a single displacement length p of three links for our 
operators. 

At this point, we can combine the smeared quark fields and gauge links into gauge- invariant 
extended three-quark operators containing the desired color and radial structure: 

®tZd*) = ^ c [?A (3) S4]fjx)[^f ) S 4 ]f,(a;)^ J Df )t 74 ]^(x) (4.8) 
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Operator type Displacement indices 



i = j = k = 
i = j = Q,k^Q 
i = 0,j = -k, k j£ 
i = 0,|j| ^ |*|,jM0 
* = -3, \j\ ^ 1*1, ^ 



Table 4.1: The spatial arrangements of the extended three-quark baryon operators 
and Quark- fields are shown by solid circles, line segments indicate gauge-covariant 

displacements, and each hollow circle indicates the location of a Levi-Civita color coupling. 
For simplicity, all displacements have the same length in an operator. 

where a specific displacement length of p = 3 has been used. An explicit factor of 74 is 
included in the creation operator ^ a p~ ^{x) to ensure that our correlation matrices are 
Hermitian. 

The displacement direction indices i,j, k = 0, ±1, ±2, ±3 define the radial structure of the 
operator. The types of displacements for the three-quark operators considered in this work 
are illustrated in Table 14. 11 In particular, comparing the effectiveness of the 
doubly-displaced-L vs. the triply-displaced-T operators can shed some light on the much 
discussed issue of Y-flux/A-flux formation of the gluon field in three-quark systems (the A 
is actually a quantum superposition of three V-flux forms). 

A sixth possible configuration in which the quarks are displaced in three orthogonal 
directions (triply-displaced-O) was not considered because the connecting gauge links for 
such an operator do not lie in the plane of the quarks. It is expected that such an operator 
would couple to the higher-lying contaminating states much more than the similar 
triply-displaced-T operator. 




Single-Site 
Singly-Displaced 
Doubly-Displaced-I 

T 

<§>-• 

Doubly-Displaced-L 
Triply-Displaced-T 
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4.3 Classification of states by transformation behavior 

ABC 

So far, we have constructed gauge- invariant extended three-quark operators <& a py,ijk( x ) an d 
®aP7,ijki x ) f rom smeared gauge links, and smeared and covariantly-displaced quark fields. 
Our goal is to use these extended operators to construct full baryon operators which 
possess good quantum numbers on our lattice. 

Quantum numbers tell us how an object transforms under some symmetry operation. As 
an example, consider the spherical harmonics Y™{Q,(fi), functions which appear frequently 
in the description of systems possessing spherical symmetry (e.g. the multi-pole expansions 
of classical electrodynamics, or the orbital wavefunction solutions of the hydrogen atom in 
quantum mechanics). The spherical harmonics Y l m (6,(p) transform under spatial rotations 
in specific /-dependent ways. The spherical harmonic F ° is a rotationally invariant function. 
The three spherical harmonics Yf 1 , Y"/ transform as a vector v having components: 

v* = ^=<YC x -Y}\ v v = -j=(Yr 1 + Yfi, v z = Yl 

These different transformation behaviors can be used to label the state of the system (in 
the example of the hydrogen atom, / is the orbital quantum number). 

We turn now to the task of combining our three-quark operators into baryon operators 
which transform irreducibly according to the symmetries of the lattice. 



4.3.1 Irreducible representations of groups 

Briefly, we define and summarize the key ideas from group theory we will need for the 
remainder of this chapter. Summation will be denoted explicitly (summation over repeated 
indices is not implied). 

A group Q is a set of elements {R} along with multiplication law which satisfies the 
following properties: 

1. Closure: RR' E Q V R, R' E Q 

2. Associativity: R 1 (R 2 R 3 ) = (RiR 2 )R3 V R 1 ,R 2 ,R-$ E Q 

3. Identity: 3 E eQ s.t. ER = RE = R V ReQ 

4. Inverse: 3 R' 1 E Q s.t. RR' 1 = R~ X R = E V R E Q 



If the multiplication operation is commutative, then the group is Abelian, otherwise, it is 
non-Abelian. A representation of a group is any set of matrices {D(R)} which satisfies the 
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group multiplication rule 

D(R)D(R') = D(RR') V R, R' G Q. 

In the following, we will denote our extended three-quark operators on each time slice r 
using a single index: 

where the master index a stands for the values of A, B, C, a, /3, 7, i,j, k, and x. 

Under a color, flavor, or spatial transformation R G Q represented by the Hilbert space 
operator U R , our operators will transform into one another as 

d 

B a -> U R B a U R = J2Bf3Wpa{R), 

(5=1 
d 

B a ^ U R B a U R = Y,BpW^{R), (4.10) 

13=1 

where d is the dimension of the (not necessarily unitary) representation matrices {W(R)}. 

For example, lattice rotations through the angles of 0,7r/2,7r, and 37r/2 about the axis Oj 
transform the different orientations of a triply-displaced-T operator into each other (taking 
Oj to be perpendicular to the plane containing the quarks). Performing the same 
operations on a single-site operator yields only one distinct operator. 

We will restrict our discussion to these types of transformations. Other types of 
transformations, such as the charge conjugation operation considered in Chapter [5], do not 
transform baryon operators into baryon operators. 

If the collection of transformations {R} being considered composes a group Q, then the 
coefficient matrices {W(R)} are a representation for Q defined on the basis {B a }. 

d 

U R ,U R B a U R Ul = J2 U R'^i3 U R' W 0a(R), 

0=1 

d d 

7=1 13=1 
d 

U RIR B a U RIR = J2B,W, a (R'R). 

7=1 

If the representation is reducible, then we can define d! < d operators 

d 

B a = ^ BpApg, a = 1, ■ ■ • ,d' < d, 

13=1 
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which transform only among themselves: 

d> 

(3=1 

where the representation matrices {W'(R)} have a smaller dimension than the original 
representation matrices {W(R)}. A representation which is not reducible is called an 
irreducible representation, or an irrep, and is labeled by a superscript in parenthesis. Such 
representations are significant because they can be used to label the different sectors of our 
Hilbert space. 

Let {D( A \R)} denote a set of riA-dimensional unitary^ matrices forming an irreducible 
representation of Q labeled by A. If we consider all nonequivalentQ irreducible unitary 
representations of a finite group Q of ordejfl g, then the quantities d\^\R) for fixed fj,, 
form a vector in a (/-dimensional space, such that: 

R V 

Eqn. 14.111 forms the backbone of our operator construction effort, and is proven in [Ham62j 
as a corollary of Schur's Lemma. The fundamental theme which will guide the operator 
construction is that of constructing operators which transform irreducibly under symmetry 
transformations. 



4.3.2 Symmetry operations 

When constructing our baryon operators, we will consider symmetry groups of the 
Hamiltonian. Consider a finite group Q of g unitary operators U R which act on the Hilbert 
space of the system (U R U R = 1). Q is a symmetry group of the Hamiltonian if each element 
commutes with the Hamiltonian: 

[U R ,H} = 0. 

Symmetry operations may be applied to operators at different times without disrupting the 
time evolution relation: 

U R B a (r 2 )U R = Une^-^B^e-^-^Ui 

= e 11 ^^ (u R B a {n)U* R ) e~ H ^\ 

We would like to find linear combinations of these operators which transform irreducibly 

3 We will be dealing with unitary representation matrices (D* (R)D(R) — 1) unless specified otherwise. 
In particular, all finite groups possess unitary representations. 

4 Two representations are equivalent if they are related by a change of basis. 
5 The order g of a finite group Q is simply the number of elements in Q . 
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under Q: 



B^ a = E c isUsB a Ul U R B^ a U R = B j>a , a Df{R) (4.12) 

s j 

Bt,a = E c i*sUsB a Ul U R B^ a U R = ]T Bf >a>a D^*(R) (4.13) 



where A labels the irrep, i labels the row of the irrep, a labels the instance of the irrep 6 ' 
and a denotes all other operator information (e.g. displacement values). 

The power of this approach may be seen by considering a general correlation function 
between two operators: 



Cgfit) = (B^ a {t)B) Ap {Q)) = Tr {B^ a (t)B), b ^) / 

= \ E Tr (4^>,a(*)4^ 6 ,/3(0)4^ 
3 R J 

= \ E Tr { U R B tMU R U R Bl b ^)U R ) (4.14) 

^ R 

= E f~ E d h ] *( r ) d V( r ) \ Tr K^K^w) 

fc,Z \ R J 

= 5Ar^-E^t,«W^(0)) (4-15) 
UA k 

where the last line came from using the orthogonality relation in Eqn. 14.111 
It is important to notice four things: 

1. The irrep and the row must be the same for the correlation function not to vanish. 

2. To get I4.14[ we invoked the cyclic property of the trace. Thus the operators are not 
orthogonal on a configuration-by-configuration basis because equation 14.151 is an 
expectation value relation. Anticipating orthogonality by calculating only 
non-vanishing correlation functions, we are effectively reducing the noise in our 
Monte Carlo estimates. 

3. The value of the correlation function is independent of the row. Thus we can get 
better statistics by correlating operators by row, and then averaging over the rows. 

4. Different instances of an irrep may mix. Because the value of the instance label does 
not affect the transformation properties of the operator, we may group the instance 
label with the auxiliary index a. Let the resulting operators be denoted B i ^ where A 
is the irrep, % is the row, and j represents the value of a and a. 



6 There may be more than one linear combination of the UjiB a U'L which transforms according to a 
particular irreducible transformation. 
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Given a set of elementary operators {B a }, we will construct the correlation matrix in a 
given irrep A as 

"-A , 



4.3.3 Projection onto irreducible representations of finite groups 

For any operators B a and B a , we construct operators which transform irreducibly under 
transformations R G Q by applying the following projection operations: 

9 R 

Bl, a = OlB a = ^J2 D ^(R)U R B a Ul (4.16) 

^ R 



The required transformation behavior is easily verified: 

U R Bly R = ^Y, D ^(S)URsB a U R s 
9 s 

= ^J2D^(R- l S)U s B a Ut 

y s 

3 



Of a acting on an B a produces an object which transforms according to the i th row of the 
irreducible representation labeled by A. Consider the effect of two such operations: 

n A n r (aw D x n (rw CUr ~b T A 



°ta°l b B a = ^J2D%>(R)Df>(S)U RS B a U 

9 R,S 

~ nA ^J2D^'(R)D$»(R-i S )U s B a Ul 

9 R,S 



k \ 1 j[ J 9 s 

= 8\r5 n n0^ h B n . 



In particular, if we define P A = 0A, we have 



P l A PfB a = 5 Ar 5 ij P l A B a 
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The P A are therefore projection operators. 



s \y a / 



= s q 

where the x (A) (^) = Tr(L> (A) ) are the characters of the group irrep, and we have used the 
character orthogonality relation |Ham62] : 



3 a 



The identity elementQ D(E) = 1 in any representation, and thus x {E) 
We can use Eqn. 14.101 to define projection matrices P^ a/3 - 



n A . 



d 



^ R 13=1 
d 



= Y, p tp B p ( 4 - 17 ) 



3 R 

In the case that the {W(R)} are not unitary, we may define the Hermitian metric matrix 
M: 



M «p=\Y,ii w u R ) w ^ R )- 



3 R 7=1 



We then compute the projection matrix for row i = 1: 

q z — ' 

y R 

From the rows of this projection matrix, r linearly-independent operators are obtained: 

d 

B l,a = E C t<*P B Pi 



7 The identity element is denoted by E, from the German einheit, meaning unit. 
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where r is the rank of the projection matrix and the superposition coefficients c^ a/3 for each 
operator are a linear combination of the rows of the projection matrix such that these 
coefficients satisfy 

d 

Yl c t<*i M i° c ip<r = <W, ol = 1, • • • , r. 

7,(7=1 

In practice, these linear combinations are obtained using the well-known Gram-Schmidt 
procedure, but with a modified inner product to incorporate the metric matrix M. The 
choice of these operators is not unique. For each of the r operators B^ a in the first row 
i — 1, we obtain the partner operators in all other rows j > 1 using: 

y 7=1 R 

4.3.4 Symmetries of the lattice Hamiltonian 

We will construct baryon operators which transform irreducibly according to the symmetry 
groups of the Hamiltonian on each time slice: 

• The SU (3) gauge group (color) 

• The Z N <S> Z N (g) Z N cyclic translation group (momentum) 

• The 577(2) isospin group (flavor) (assumes m u = m d ) 

• The 0® crystal point group (spin and parity) 

4.3.5 Color 

By requiring our operators to be gauge-invariant, we have already ensured that the 
operators transform irreducibly as color singlets. The irreducible representation of the color 
group by which the operators transform is the trivial representation where 

D(R) = 1 V Re SU(3). 

An the projection operation for this representation acts on the color indices only, and 
results in the Levi-Civita coupling e a bc- 
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4.3.6 Momentum 



The next index we consider for projection is the spatial reference site x. Spatial 
translations form an Abelian group (the product of two translations is a translation, each 
translation possesses an inverse, et cetera). The displacement operator acting on Euclidean 
space-time functions is given by 



ia-iP. 

3 J 

D(a)f(x) = f(x-a). 



D(a) = e m ^\ Pj = idj 



Here D(a) is an operator, not a representation matrix. However, this abuse of notation will 
not lead to trouble because the irreps D^ k \a) of the group of translations are all 
one-dimensional and given by: 

D {k \a) = 

where kj is the irrep label (the momentum). 

Functions which transform irreducibly under translations are labeled by their momentum k: 

f( k )( x ) = ce~ ik ^ 
D(a)f {k) (x) = f {k \x -a) = f {k \x)D {k \a) (no sum on k) 
PJ (k) (x) = kj^(x) 

For our operators, we have a finite number Nj of discrete spatial lattice sites in each 
direction j with periodic boundary conditions: 

D {k \N jas ) = D {k) (0) 

where a s is the spatial lattice spacing. Thus, our group of spatial translations is a direct 
product of three cyclic groups Z N . = {z , z 1 , - • • , z Nj }, where 

• z is the identity element 

• z\ corresponds to a translation from x to x + j 



N 3 

The unique allowed momenta on the lattice are given by the irrep labels of Z^r. 

kj = — — rij, (no sum on j), rij — 0, 1, • • • , Nj — 1, . 

jOjg 

Some researchers shift the range of rij to allow negative momentum values. Because we are 
interested in the mass spectrum, we will work in the zero-momentum sector 
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k\ = k2 = k 3 = 0. Applying the projection operation (Eqn. I4.16P with the trivial 
representation gives the average over all spatial sites. Therefore, our zero-momentum 
extended three-quark operators are translationally invariant and given by: 

= jjj^- £ e*® D i r)®Df%(x, t)®D®%(x, r) 74 (4.19) 

X 



4.3.7 Flavor 



For this work we will work with 'light baryons,' containing u, d, and s valence quarks only. 
Correspondingly, our flavor indices A, B, and C, can only take on the values u, d, s. In 
addition to conserving flavor, our Hamiltonian possesses an SU (2) isospin symmetry in 
which the u and d quarks are assigned the same mass (m u = m^) and are treated as 
components of an isospin doublet: 

/ u 



I 



d 



where / stands for light quark. The group of transformations on the isospin doublet is 
SU(2), the continuous group of all 2 x 2 unitary matrices with determinant = +1. 

In nature isospin symmetry is broken by a small u-d mass splitting which arises due to 
electromagnetic interactions and other effects, but this splitting is less than 1% mass of the 
lightest baryon, the proton. It is therefore reasonable to treat SU(2) isospin symmetry as 
exact in lattice QCD simulations at present levels of precisions. We will require our 
operators to transform irreducibly under isospin transformations. 

On the other hand, the action explicitly breaks SU (3) flavor symmetry (m u = = m s ) 
by the designation of a heavier strange quark mass. In nature, the strange quark - light 
quark splitting is about 15 — 30% of the proton mass. Although SU(3) flavor symmetry is 
not exact, we may design our operators to transform according to SU(3) flavor irreps in 
addition to SU (3) isospin irreps, but must expect such operators to mix. 

The underlying mechanism causing the different SU(3) irreps to mix is illustrated by 
considering the effect of using the projection operation (Eqn. 14.161) for a group which is not 
a symmetry of the Hamiltonian: 

~D A ' a ^\ _ Ua n( A )* 



BZ(r) = -^J2D^(S)U s B a (r)Ul 



9 s 



J E D^(S)U s e +H ^B a (0)e-^Ut 



s 



If [U R , H) 0, then 
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and the orthogonality relation 14.151 does not hold. If we require our baryon operators to 
transform according to the S77(3) flavor group, then the explicit symmetry breaking 
(m s 7^ m u ,m^) in the Hamiltonian will cause the operators in different symmetry sectors to 
mix. This mixing will occur among baryon operators with equal isospin I, isospin 
projection ^3, and quark flavor content (e.g. strangeness). Some different baryon sectors 
labeled by the good quantum numbers I and S are illustrated in Table 14.21 



I 


s 


Baryon Sector 


1 

2 





N°,N + 


3 
2 





A" A , A+, A++ 





-1 


A 


1 


-1 




1 

2 


-2 







-3 





Table 4.2: The isospin and strangeness quantum numbers for the different baryon sectors. 
The isospin projection I3 ranges from —I to I by increments of 1. The charge Q of each 
baryon listed increases with increasing isospin projection J 3 . 

SU (2) isospin is a subgroup of SU (3) flavor, and consists of the set of SU(3) flavor 
elements which leave the s quarks invariant. Given an irreducible representation of SU(3) 
flavor (which is not necessarily three-dimensional), the subset of irrep matrices 
corresponding to the SU{2) subgroup elements forms a representation of the SU(2) isospin 
group. This representation will generally be reducible, which implies that we can find a 
basis in the SU(3) irrep such that the SU(2) representation (which is not two-dimensional) 
can be reduced to block diagonal form. In this basis, we will have states which transform 
irreducibly under both SU(3) flavor and 577(2) isospin. We will return to this theme in 
section 14.4.41 when we apply this method of subduced representations by restricting the 
continuum spin group 0(3) to the finite lattice rotation group. 

4.3.8 Flavor irreps 

The finite group projection operation (Eqn. I4.16P is readily extended to compact 
continuous groups^ such as SU(2), but it is simpler to use Lie algebra methods to construct 
operators which transform irreducibly. We will write out sums explicitly in this section 
(summation over repeated indices is not implied). 

8 Elcments of a continuous group, such as the group of three dimensional rotations, are characterized by 
a set of real parameters. If there exists a closed set which covers the complete range of the parameters, 
then the group is said to be compact. Elements of U(l) may be written as exp(i0), where the parameter 
9 G [0, 2n) C [0,2tt]. Thus U(l) is a compact continuous group. 



87 



An element R of SU (2) can be written as 

3 

R = exp(i a i ( R ) T j), t j = y 

where the Oj are the traceless Hermitian Pauli spin matrices as defined in Eqn 12.101 The 
SU (2) group structure is determined by the commutation relations among the generators: 

[n,Tj] = i ejjkTk 
k 

We can construct a Casimir operator which commutes with all of the group generators: 

r 2 = T 2 + T 2 + T 3 2 , [r 2 ,Ti] = (i = l,2,3). 

By Schur's Lemma, we know that r 2 is proportional to the identity matrix and can thus be 
used to label the different irreps. We may also choose the basis of Hilbert space states in 
each irrep such that r 3 is diagonal. By introducing the ladder operators 

r± = 7i ± ir 2 

which satisfy 

fa, t±] = ±r±, r 2 + -r+r_ + -r_r + = r 2 , 
it can be shown (c.f. jMer98j ) that 

r±\I,h) = ^ I {! + !)- h{h± 1 )| I, h±l) 
= y/{lThW±h + l)\I,h±l) 

T 3 \I,I 3 ) = I 3 \I,I 3 ) 
T 2 \I,I 3 ) = I(I+1)\I,I 3 ) 

where the different irreps are labeled by I and the row is labeled by I 3 . 

By definition, the irreducible representation matrices {D^} (also known as Wigner 
rotation matrices in this context) are the matrix elements of the SU(2) operators between 
different irrep basis states in the Hilbert space: 

u R \i,h) = Y^\iJ's) D ^ h (R) 

D% h (R) = (I,%\U R \I,I S ). 

We would like our baryon operators to transform irreducibly under SU(2) isospin 
t r ansf or mat ions : 

= ^^,(j,4i^|/,j 3 ) 
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where we have suppressed all indices except the flavor irrep and row. Considering an 
infinitesimal transformation dj — > (i = 1, 2, 3) and keeping terms up to first order in a, at 
each step gives: 

j k I' 3 I 

3 j I' z 



where in the last step we used the fact that the infinitesimal parameters are arbitrary. 
Our baryon operators will transform under isospin according to the irreducible 
representation labeled / if and only if 

fo-Xl = V(/-/8)(I + /S + I)5j, +X) 



From these relations, we see that the analog of the Casimir operator is: 

fa, [r 3 ,K]] + \[r + , [r-,K]] + \\r-, [r+.sj,]] = /(I + l)Sf a . 

We may now proceed to construct baryon operators from the individual quark fields. The 

"it _ ~ d 

light quarks are represented by u = ?/> and d = ijj smeared source operators having 7=1 

~ s 

and ^3 = +1/2, —1/2 respectively. The strange quark is represented the s = ip smeared 
source operator having I = 0, 1 3 = 0. These operators therefore satisfy 

[t 3 ,u] = \u, [r 3 ,J] = [r 3 ,s] = 0, 

[t+,U] = 0, [r + ,d] = m, [r + ,s] = 0, 
[t_,u] = d, [t-,^] = 0, [t~-,s] = 0. 

Due to isospin symmetry in the Hamiltonian, our particle masses will not depend on the 
irrep row I 3 . We will therefore construct only one operator in each isospin /, strangeness S 
sector, the operator corresponding to maximal isospin projection ^3 = /. We now have 
three- quark elemental operators having definite isospin /, maximal I3 = I, and strangeness 
S: 

A,B,C 

where ^\bc are the coefficients for the flavor channel F. The three-quark elemental baryon 
operators are given explicitly in terms of the extended three-quark operators in Table 14.31 
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Baryon I — h S Operators 

' A++ H Q 

y+ i _i ^" MS 

^ 11 ^a/3-y,ijk 

^ -^uud -rr-duu 

2 ^a/3y,ijk ^a/3^/,ijk 

-o i o $ ss " 

1 — 1 2 a/3j,ijk 

AO n i ^ urfs ^T^" 5 

— F 

Table 4.3: Elemental three-quark baryon operators B a/3 ^ having definite isospin J, maxi- 
mal J 3 = 7, and strangeness S in terms of the gauge-invariant extended three quark operators 
7f> ABC ( \ 



4.4 Spin and parity 



At this point, we have used smeared quark fields and gauge links to construct three-quark 
elemental operators B i and Bf , which have the desired momentum, color, flavor, and 
radial structure. The superscript F denotes the flavor structure, and can take on the values 
(A), (E), (TV), (S), (A), or (fl). The index i now consists of only the Dirac spin indices 
a/3'y and the displacement direction indices ijk. Our final task in constructing our 
operators is to project our operators onto states of definite spin and parity. 



4.4.1 Baryon spin in the continuum 

In the continuum, we could construct our baryon operators using the irreducible 
representations of the continuous spin group 577(2), the spinorial double cover of the 
rotation group 573(3). The irreps of SU(2) are labeled by the spin quantum number 

j _ 1 3 5 

2 ' 2 ' 2 ' ' ' ' ' 

For example, we can construct a single-site spin-1/2 nucleon baryon operator by taking 

bL n \t) = j:b^ j ^ (t)^ 

Pi 

where T is some 4x4 matrix of constants. 

By combining two of the quarks into a Lorentz scalar or pseudoscalar, the resulting 
three-quark operator transforms like a spin-1/2 Dirac field under Lorentz transformations. 
If two of the quarks are combined into a four- vector (or axial four vector), then the 
resulting interpolating three-quark operators produce both spin-3/2 and spin-1/2 states. 

Although this method of building up baryon operators is in widespread use, it becomes 
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extremely cumbersome when constructing higher spin states or complicated extended 
operators. The difficulty is compounded by the fact that the lattice regulator explicitly 
breaks Lorentz covariance. 



4.4.2 Baryon spin on the lattice 

The discretization of space onto a lattice breaks the rotational symmetry of the system. 
Thus, different continuum spin states will mix, making the identification of those states 
problematic. We may mitigate this difficulty, however, by designing our operators to 
transform under the lattice symmetry group Of, the double- valued (spinorial) octahedral 
crystallographic point group. The states excited by our operators will have well-defined 
lattice spin and parity labels, and will not mix with each other. 

Applying the projection operation (Eqn. 14. 16j) gives the full baryon operators: 

\ T ) — ^af}f,jkl C a3^,jkV (4 2T) 

tjAXF ( \ — r>F _(AAiJ* V*-* 1 -) 

where A is the Of irrep and A is the row. The coefficients c(AAz) aj g 7 jw can be determined 
using the methods of Subsection 14.3.31 once we have the irreducible representation matrices 
for Of. 

The exposition below is adapted from [BEF + 05b] and [MEF + 04] . 



4.4.3 Irreducible representation matrices for Oj^ 

The basic building blocks used to assemble our baryon operators transform under the 
allowed spatial rotations and reflections of the point group Of according to 

uM{x)DfA A Ui = (1p{Rx)D^) A S Pa (R), (4.22) 

Unffifa))^ = S{R)-l(b%i,{Rx)) A ap , (4.23) 

where the transformation matrices for spatial inversion I s and proper rotations C n j through 
angle 2Tr/n about axis Oj are given by 

S(C nj ) = exp(§av[7^7j), (4-24) 
S(I S ) = 74, (4.25) 

with LUki = —ZTfEjki/n and oj^ = uiki = [oj^ is an antisymmetric tensor which 
parameterizes rotations and boosts). A rotation by 7r/2 about the y-axis is conventionally 
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denoted by C^ yi and C^ z denotes a rotation by 7r/2 about the z-axis. These particular 
group elements are given by 

S(C % ) = -L(l + 7l73 ), 5(C 4 ,) = -^(l+ 7 27i)- (4-26) 

The allowed rotations on a three-dimensional spatially-isotropic cubic lattice form the 
octahedral group O which has 24 elements. Inclusion of spatial inversion yields the point 
group Oh which has 48 elements occurring in ten conjugacy classes. All elements of Oh can 
be generated from appropriate products of only C^, C^, and I s . 

Operators which transform according to the irreducible representations of Oh can then be 
constructed using the group-theoretical projections given in Eqn. (14.211) . Orthogonality 
relations and hence, projection techniques, in group theory apply only to single-valued 
irreducible representations. However, the fermionic representations are double-valued 
representations of Oh- The commonly-used trick to circumvent this difficulty is to exploit 
the equivalence of the double- valued irreps of Oh with the extra single- valued irreps of the 
so-called double point group 0®. This group is formed by introducing a new element E 

2 

which represents a rotation by an angle 2n about any axis, such that E = E (the 
identity). By including such an element, the total number of elements in Oj^ is double the 
number of elements in Oh- The 96 elements of 0® occur in sixteen conjugacy classes. 

Since baryons are fermions, we need only be concerned with the six double-valued irreps of 
Oh- There are four two-dimensional irreps Gi g , Gi u , G 2g , and G 2u , and two 
four- dimensional irreps H g and H u . The subscript g refers to even-parity states, whereas 
the subscript u refers to odd-parity states^. The irreps G\ g and G\ u contain the spin-1/2 
states, spin-3/2 states reside in the H g and H u , and two of the spin projections of the 
spin-5/2 states occur in the G 2g and G 2u irreps, while the remaining four projections reside 
in the H g and H u irreps. The spin content of each Oh irrep in the continuum limit is 
summarized in Table 14.61 This table lists the number of times that each of the Oh irreps 
occurs in the J = |, §,§,•' ' representations of SU{2) subduced to Oh- 

To carry out the projections in Eqn. (14.211) . explicit representation matrices are needed. 
Our choice of representation matrices is summarized in Table 14.41 Matrices for only the 
group elements C^C^, and I s are given in Table |4~41 since the representation matrices for 
all other group elements can be obtained by suitable multiplications of the matrices for the 
three generating elements. For baryons, the representation matrix for E in each of the 0® 
extra irreps is —1 times the identity matrix. 

To give an example of this method, the single-site N + (nucleon) operators which transform 
irreducibly under the symmetry group of the spatial lattice are given in Table 14.51 



From the German gerade meaning 'even,' and ungerade meaning 'odd' 



92 



A 



r( A )(c, 



'hj> 



T {K \G 



lz I 



G 



G 







1 


' 1 


-1 








1 


1 






-1 


" 1 


-1 






7! 


1 


1 




1 




-1 


1 




-1 - 


1 






Vs 


1 - 


1 - 


-Vs 




1 


V3 


1 



1 

71 



1 

71. 

-1 " 
71. 



l-z 


l-z 




-1-i 
1-i 
l+i 
-1 




l+i 



l+i 







Table 4.4: Our choice of the representation matrices for the double-valued irreps of Oh- The 
Giu, G2U1 H u matrices for the rotations C^C^ are the same as the Gi g , H g matrices, 
respectively, given below. Each of the Gi g , Gig-, H g matrices for spatial inversion I s is the 
identity matrix, whereas each of the Gi u , Giu-, H u matrices for I s is —1 times the identity 
matrix. The matrices for all other group elements can be obtained from appropriate multi- 
plications of the C 4y , C43, and I s matrices. 



Irrep 


Row 


Operator 


Row 


Operator 


Gi g 


1 


#211 


2 


#221 


Gig 


1 


#413 


2 


#423 


Gig 


1 


2^332 + ^413-2^431 


2 


2^432-2^444-^423 


Giu 


1 


#433 


2 


#443 


Giu 


1 


#321 — #312 


2 


#421 — #412 


Giu 


1 


#312 + #321-2AT 411 


2 


2iV322 — #412 — #421 


Hg 


1 


^#331 


2 


#332-#413 + 2#31 


Hg 


3 


2^432 + ^444-^423 


1 


\/3#442 


H u 


1 




2 


—#312 — #321 — #411 


H u 


3 


—#322 — #412 — #421 


4 


— #422 



Table 4.5: The single-site N + operators which transform irreducibly under the symmetry 
group of the spatial lattice, defining N af 3 y = $q^ t>000 — $a/3r,ooo ( see Table (14. 3p ) and using 
the Dirac-Pauli representation for the Dirac gamma matrices. 

4.4.4 Continuum spin identification 



Once we have the lattice spectrum, we may use the method of subduced representations to 
to identify the continuum spin values of the extracted states. Because is a subgroup of 
the continuum rotation group for spinors (with inversion), the continuum states appear as 
degenerate states within each of the lattice irreps. Table 14.61 denotes the number of times a 
continuum-limit spin J state will appear in each irreducible representation of the lattice. 
Thus, Table 14.61 allows us to identify the continuum-limit spin J corresponding to the 
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masses extracted in our Monte Carlo calculations. For example, to identify an even-parity 
baryon as having J = 1/2, a level must be observed in the G\ g channel, and there must be 
no degenerate partners in either of the G2 g or H g . A level observed in the H g channel with 
no degenerate partners in the G\ g and G23 channels (in the continuum limit) is a J = 3/2 
state. Degenerate partners observed in the G2 g and H g channels with no partner in the G\ g 
channel indicates a J = 5/2 baryon. In other words, Table H~6l details the patterns of 
continuum-limit degeneracies corresponding to each half-integral J value. 



J 






n k 


1 

2 


1 








3 
2 





1 





5 
2 





1 


1 


7 
2 


1 


1 


1 


9 
2 


1 


2 





11 
2 


1 


2 


1 


13 

2 


1 


2 


2 


15 
2 


1 


3 


1 


17 

2 


2 


3 


1 



Table 4.6: Because is a subgroup of the continuum spin group, the continuum states 
appear as degenerate states within each of the lattice irreps. This table allows us to identify 
the continuum spin state which corresponds to each lattice state. 



Chapter 5 

Evaluation of baryon correlation 
matrices 

5.1 Charge conjugation and backward propagating 
states 

5.1.1 Charge conjugation 



As discussed in chapter |5J the Euclidean ip and ip operators are the analytic continuations 
of ip M and ipM, respectively, and must be treated as independent fields: 

if) ^ V f 74- 

However, we can define a charge conjugation operator Uq on the Hilbert space which 
exchanges quarks and antiquarks. This introduces a relation between the operators ip an d 



where the 4x4 charge conjugation matrices C and C keep track of how the spin 
components of the fields change under charge conjugation and summation over repeated 
indices is implied. To find the form of C and C, we consider how the fermion terms in the 
continuum Euclidean action transform under the action of charge conjugation: 



Uc4> a (x)W c = ippCp, 
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where we anticommuted the Grassmann variables, integrated by parts, and discarded the 
surface term. In order for the action to be invariant under charge conjugation, we must 
have 

cc T = -1, C 1 »C T = 7 J, 

A,{x) -> -Al(x), U,{x) -> Ul(x). 

A common choice for the charge conjugation matrix in the Dirac-Pauli basis (c.f. 
Chapter M) is: 

= C = 7472, 

and the fields transform as: 

i) -> ^7472, ^ -> ^7472- 

As an example, consider the charge conjugation behavior of the simple bilinear form ip^^ip 
(which can be used to interpolate for vector mesons): 

i>li$ -> ^74727m(7472) T ^, 
= ^74727^7274^, 
= -^7^, 
= i>l^- 

5.1.2 Backward propagating states 

Our operators are constructed using fermion fields ip(x) which create a quark and 
annihilate an antiquark. Hence, each of our baryon operators creates a three-quark system 
of a given parity P and annihilates a three-antiquark system of the same parity P. This 
means that in the baryon propagator, a baryon of parity P propagates forward in time 
while an antibaryon of parity P propagates backwards in time. If we apply charge 
conjugation to our baryon operators, we arrive at the 'backward propagating' partners 
mentioned in Chapter [2j 

Unlike boson fields, a fermion and its antifermion have opposite intrinsic parities, so that 
the antibaryon propagating backwards in time is the antiparticle of the parity partner of 
the baryon propagating forwards in time. Since chiral symmetry (vanishing quark mass) is 
spontaneously broken in QCD, the masses of opposite-parity partners may differ. The 
forward propagating baryon will have a mass different from that of the antibaryon 
propagating backwards in time. If the even- and odd-parity baryon operators are carefully 
designed with respect to one another, it is possible to arrange a definite relationship 
between the correlation matrix elements of one parity for r > and the opposite-parity 
matrix elements for r < 0, allowing us to increase our statistics. 
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We can get a useful relation between the correlation matrix elements C ab (r) and C ab (T — r) 
if we utilize the projection procedure in Chapter @] for our ei>en-parity operators, and utilize 
charge conjugation to construct our odd-parity operators. Consider the correlation matrix 
element of two even-parity operators for r > 0. Suppressing irrep, displacement, and flavor 
indices, one sees that invariance under charge conjugation implies that 

= c^ 7 cg_(0| C^CB a ^rpCB^(0pC |0>, 
= cS; C g_<0| S QW (r)S^ v (0) |0> 

= c«; C g_(0| S QW (0)S^ v (-r) |0) 

2 2 2 2 2 2 

using invariance under time translations of the above expectation value and invariance of 
the vacuum under charge conjugation. The last line above represents the correlation of 
odd-parity operators propagating temporally backwards. Hence, for a given even-parity 
operator B?(t), we can define an odd-parity operator B^(t) by rotating the three Dirac 
indices using the 72 matrix and replacing the expansion coefficients by their complex 
conjugates such that the correlation matrices of the even- and odd-parity operators are 
related by 

C G ab l3 {r) = -CS ta (-r)*, 

C G ab 2g {r) = -Cf-(-r)*. 

For a lattice of time-extent T with anti-periodic temporal boundary conditions, this means 
that 

C { f b \r)=Cl:\T-r), (5.1) 

and similarly for the other irreps. This allows us to appropriately average over forward and 
backward temporal propagations for increased statistics. 

5.1.3 Improved estimators 

In the absence of any external applied fields, the energies of the baryons do not depend on 
the row A of a given irrep A, so we can increase statistics by averaging over all rows: 




Because ip and if) are independent fields in the functional integral, the correlation matrix 
for a single configuration will not be Hermitian due to statistical fluctuations inherent in 
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the Monte Carlo estimate. Thus, we can increase our statistics by taking the average: 



1 



*(<?), 



In addition, we can use the backward propagating states in our correlation matrix to 
improve the statistics of the opposite-parity correlation matrix by exploiting the relation 



cS(r) = CF a >(T-T 



(«), 



Putting it together gives us the improved estimator for our correlation matrix elements: 



1 

4 
1 



CS(T) + C£'(T - r) + C™» + C^*(T - r) 



i{u). 



-1(9)* 



'(«)*/ 



C^{r) + CZ\T - r) + C^(r) + C%*(T - r 



f(s)*/ 



(5.2) 
(5.3) 



It is important to note that although we have increased the precision of our estimates by 
averaging, we have not increased the number of configurations produced by the Monte 
Carlo method. Thus, when estimating the jackknife errors according to the method 
presented in Chapter [3j we must form our n th jackknife sample before calculating the 
improved estimators. Specifically, we do not treat the different terms in the improved 
estimator (Eqn. 15. 3j) as additional configurations which can be excluded to form a jackknife 
sample. 



5.2 Three-quark propagators 



Having constructed our baryon operators, we now have correlation matrix elements which 
look like: 




1 n A Mb 

= 7rEE c r c r^(^(o)). 

A A=l k,l=l 

where the indices a, b, k, I range from 1, ■ • • , Mb, and encapsulate all of the information 
about the type of operator (single-site, singly-displaced, doubly-displaced-L, 
doubly-displaced-I, or triply-displaced-T), and the particular linear combination over Dirac 
indices used to make an operator which transforms irreducibly according to row A of the A 
irrep (Gi g , H g , G 2g , G\ u , H u , or G 2u ) in the flavor channel F (A, E, N, S, A, or fl). 

Since the number of elemental operators is large and the quark propagators are rather 
expensive to compute, it is very important to use symmetry to reduce the number of 
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Operator Type 



Displaced Quark Direction(s) 



Singly-Displaced 



+z 
+z, -z 

+y, +z 

+z, +y, -z 



Doubly-Displaced-I 

Doubly-Displaced-L 

Triply-Displaced-T 



Table 5.1: To minimize the number of sources, thereby reducing the number of three-quark 
propagators needed, we rotate our source and sink operators such that the displaced quarks 
at the source are always in the same canonical positions. 

quark-propagator sources. Given the cyclic property of the trace and the unitarity of the 
symmetry transformation operators, we know that 



for any group element R of Oh- Hence, for each source B l (0), we can choose a group 
element Ri such that we minimize the total number of source elemental operators which 
must be considered. For example, consider the singly-displaced operators. We can choose 
an Ri such that the displaced quark in the source is always displaced in the +z direction. 
Similarly, a group element Ri can always be chosen to rotate each of the other types of 
operators into a specific orientation. The canonical source orientations for each type of 
operator are given in Table 15.11 

To make zero momentum operators, we averaged over all spatial reference sites x at both 
the sink and the source. However, we know that correlation functions of operators having 
different momenta will vanish (as expectation values), and thus we need only perform the 
spatial average at the sink. The source operator, evaluated at only one reference site, will 
no longer be translationally invariant, but the zero momentum contribution will be the 
only one which gives a non-canceling contribution to the correlation function. This 
increases the noise in our Monte Carlo estimates, but dramatically reduces the number of 
quark propagators needed. 

The coefficients cf^ F in the baryon operators involve only the, Dirac spin components and 
the quark displacement directions and are independent of the color indices and spatial 
sites. Thus, in calculating the baryon correlation functions, it is convenient to first 
calculate gauge-invariant three-quark propagators in which all summations over color 
indices and spatial sites have been done. A three-quark propagator is defined by 



(BF(t)B? (0)> 



(U R B^(r)U R U R Bf (0)4), 



^W Vk (R)W V i(Ry(B^T)Bf,{0)) } 



k',l'=l 



q(ABC_)(pp)_ 

(ai|m)(/3j|/3j)(7fc|7fc) 



(*) 




x 



X 




(5.4) 
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where G^.^_-r_(x,t\xo,0) denotes the propagator for a single smeared quark field of flavor 
A from source site xq at time t — to sink site x at time t. At the sink, a denotes color, a 
is the Dirac spin index, % is the displacement direction, and p is the displacement length, 
and similarly at the source for a, a, i, and p, respectively. Notice that the three-quark 
propagator is symmetric under interchange of all indices associated with the same flavor. 
As usual, translation invariance is invoked at the source so that summation over spatial 
sites is done only at the sink. These three-quark propagators are computed for all possible 
values of the six Dirac spin indices. 

Each baryon correlation function is simply a linear superposition of elements of the 
three-quark propagators. These superposition coefficients are calculated as follows: first, 
the baryon operators at the source and sink are expressed in terms of the elemental 
operators; next, Wick's theorem is applied to express the correlation function as a large 
sum of three-quark propagator components; finally, symmetry operations are applied to 
minimize the number of source orientations, and the results are averaged over the rows of 
the representations. C++ code was written by Colin Morningstar to perform these 
computations, and the resulting superposition coefficients are stored in computer files 
which are subsequently used as input to the Monte Carlo runs. 

Wick's theorem is an important part of expressing the baryon correlation functions in 
terms of the three-quark propagators. To simplify the notation in the following, let the 
indices fi, u, r each represent a Dirac spin index and a displacement direction, and suppress 
the displacement lengths. Define c^ T = c^*, T ,7^,7^,7^ r ,, then the elements of the baryon 
correlation matrix in the A ++ channel are given in terms of three-quark propagator 
components (before source-minimizing rotations) by 

r^(t\ - / r { - uuu) (t\ 

U ij \ L ) ~ C ^r C Tn^\ U (r\jl)(u\V)( fJ ,\T)y i ) 

^(uuu) ,,x ^(uuu) /,x 

"T U {r\-p){ V \T)^\V)\ l ) "I- ^(T\v)(v\-p)(n\T)\ l ) 
^(uuu) /,x ^{uuu) /,x 

+ U (t|F)Mt)( M |7I)W ^ u (t|-)Hf)(Hm)^ 

The N + correlation functions are expressed in terms of components of three-quark 
propagators by 



(uud) 

Oi|jl)(l/|P)(T|T) 



+ ^(p\u){u\Jl){r\-) ~ ^jn\f){u\V)(r\Tl) 
^(Hf)H-)(t|m) ~~ ^MF)(t|JZ)(m|t) 
^(u\JI)(t\V)^\t) + U (t\t)(v\V)(h\TZ) 
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and for the S + and A channels, one finds 

C^ff) - r {i) r (j) fa {uus) (f) 



+ G Kmp)(t|t)(*)}> ( 5 - 7 ) 

r < - A U+\ - 7:0') (P<(uds) 

^ij i l ) - C ^T C Jli7fy- r (p,\Tl)(u\V)(T\T) 



f-i(uds) fi{uds) 
U (»\vXv\JZ)(t\t) ~ ^(u\Tl)(^)(r\r) 



^ti)MMr\r)}- ( 5 - 8 ) 



(0\Tf(^U) |0> = 



5.2.1 Smeared, covariantly-displaced quark propagators 

We are now ready to compute the quark propagators for smeared and displaced sinks and 
sources. Recall that we are computing expectation values via path integrals, such as 

J V@, V, U) M V, U) exp (-mUty - S G [U\) 

J V&,4>,U) ex P (-^Q[U}^-S G [U]) 

VU h($,ij),U) detQ[U] exp(-S G [U}) 

J VU det Q[U] exp(-S G [U]) 

where U represents the gluon field (link variables) and <Sg[C/] is the pure gauge action. The 
remaining path integration over the link variables is carried out by the Monte Carlo 
method. Results in the so-called quenched approximation are obtained if det<5[£/] is ignored 
in the updating and measurement process. Full QCD simulations require the incorporation 
of the computationally-expensive fermion determinant detQ[C/] in the updating. 

Now define 



Vty,i>) (Df^) Aaa (x,t) ^Df)^(x O: 0) exp {-^Q[U]^) 

{A) , 
aaip\aaip^ 



= detQ[£/] 5 A - A G^Jx,t\x ,0). 



To see how to compute the quark propagator G^.,_j_(x,t\xo,0), let us proceed in steps of 
increasing complexity. First, we shall focus on the application of the three-dimensional 
Laplacian. If we can understand how to compute the quark propagators with the 
three-dimensional covariant Laplacian acting on both the source and sink fields, then it is 
straightforward to understand how to compute the smeared-smeared propagator. Finally, 
the p-link displacements will be included. 



101 



First, consider the application of a single Laplacian on both the source and sink fields: 

(AV)» = E {v?(x)Mx+3)-Mx)}, (5-9) 

j=±l,±2,±3 ^ J 

(M) (*) = E (^(x+j)C/j ba (x)-^ a (x)}, (5.10) 

j=±l,±2,±3 ^ J 



where the indices a, 6 are color indices. To simplify the notation, both flavor and spin 
indices are omitted, as well as the gauge-field smearing tildes. Our goal now is to compute 

J PtfM) (A^) a (x) (^K) b (y) exp Q[C/] V) • (5.11) 

First, write the Laplacian operators in matrix form: 

A ab (x\y) = £ [c/f(x)5(x+3-,y)-^(x,y)}, 

= Ej^ ^ <^+^) + U} a \x-j) S(x-j : y) - 25 ab 5(x,y) 

3=1 ^ 

= E^f^) 5(x+j,y) + C/f(y)* *(x-3-,y)-2^(x,y)j, (5.12) 
A a6 (x|y) = E \uj ab (y)S(x,y+j)-5 ab 5(x,y)\, 

= ll\ U } ab (y) + Uf(y-j) S(x,y-j) - 25 ab 5(x,y) 



3=1 K 

= Y,\ U 3 a (yy S(x-ly) + Uf(x) 5(x+j,y)-25 ab 5(x,y)\, 

3=1 ^ ' 

= K b (x\y). (5.13) 
Grassmann integration quickly gives us the quark propagator: 

J V) (AV) a (x) (fflfi) exp (-^ Q[t/] V) (5.14) 

= y P(^, V) A ac (x|z) ^ c (z) A d6 (^|y) exp (-^ Q[C/] (5.15) 

= detQ[E/] A ac (x|z) Q-J{z\w) A db (w\y). (5.16) 

A straightforward application of the above formula would require the computation of the 
inverse of the Dirac matrix Q for several local sources d. For higher powers of the 
Laplacian which may arise in usual quark smearings, the number of values of d needed may 
be quite large. Due to the computational expense of this inversion, it is better to proceed 
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as follows. The inverse of the fermion Dirac matrix is computed by solving the linear 
system of equations 

Qf c {v\z) Q~ d {z\w) = S fd 8(v,w). (5.17) 
Apply A from the right onto both sides of this equation to obtain: 

Q fc (v\z) Q^(z\w) A db (w\y) = 6f d 5(v,w) A db (w\y). (5.18) 

Write 

P cb (z\y) = Q~J{z\w) A db (w\y) 7 (5.19) 
R fb (v\y) = 5 fd 5(v,w) A db (w\y), (5.20) 

then the above equation becomes 

Q fc (v\z) P cb {z\y) = R fb {v\y). (5.21) 

This is a linear system of equations. We can solve for P cb (z\y) for a single fixed value of b 
and y by some variant of the conjugate gradient method. This is then repeated for the 
different colors b and Dirac spin components a, but usually only one site y is involved. The 
final propagator is obtained using 

G ab {x\y) = A ac {x\z) P cb {z\y). (5.22) 

Eq. (15.20p tells us how we must apply the Laplacian for the source, and Eq. ( 15. 22ft tells us 
how we must apply the Laplacian for the sink. 

The generalization of these results to higher powers of the Laplacian is straightforward. 
One starts at the source by forming the right-hand side 

Rab(x\y) = 5 aci 8(x, zi) A ClC2 (zi\z 2 ) (z 2 \z 3 ) ■ ■ ■ A Cnb (z n \y), 

= A aCn (x\z n ) . . . A C3C2 (z 3 \z 2 ) A C2Cl (z2\z!) 8 cib 5(zi,y), (5.23) 

solves the linear system of equations 

Q ac {x\z) P cb {z\y) = R ab {x\y), (5.24) 

then applies the Laplacians at the sink to the result: 

G ab {x\y) = A )... A C:iC2 (z 3 \z 2 ) A C2Cl (z 2 \zi) P cib {zi\y). (5.25) 

Using Eq. (15. 12j) . the application of the Laplacians at the sink in Eq. (15.251) is 
straightforward. For fixed b,y, P cib (z\\y) can be viewed as a known vector (since it is a 
column of a matrix), then the applications of the Laplacians are equivalent to successive 
matrix-vector multiplications. Similarly, using Eq. ( 15.131) in Eq. (I5.23f) at the source is 
straightforward, too. 

The computation of the propagator for smeared fields 

J V$,i>) 4> a (x)^b(y) exp (-0 Q[U] V) (5.26) 



103 



is now obvious. Let the smearing kernel be 

F ab (x\y) = 5 ab 5(x,y)+£ A ab (x\y), (5.27) 

= (l-6Z)5 ab 5(x,y)+Zj2{u?(x) 5(x+j,y) + U^(y)* 5(x-j, y) 1.(5.28) 

j= i < J 

One first forms the source field using 

R ab (x\y) = F%(x\z) 5 cb 5(z,y), (5.29) 
= F aCni (x\z ni ) . . . F C3C2 (z 3 \z 2 ) F C2Cl (z 2 \zi) 8 Clh 8(z u y), (5.30) 

solves the linear system of equations 

Qac(x\z) Pdb{z\y) = Rab(x\y), (5.31) 

then applies the smearing at the sink to the result: 

G ab {x\y) = F£(x\z) P cb (z\y), (5.32) 

(z 3 \z 2 ) {z 2 \zi) P Clb (zi\y). (5.33) 

Lastly, we need to include the p-link covariant displacement operators: 

J V@, (p¥>$) a (x) (Vrfp) ^ (y) exp (-^ Q[U] ^) • (5.34) 

In matrix form, these operators are given by 

(D^Jx^) = Up(x)U^(x+j)...U; p \x + (p-l)j)5(x+pj,y), (5.35) 

W) iAv) = ^r p (y+(p-i)jO---^ C2 (y+jO^ C26 (i/)5(^i/+pj), (5.36) 



j 

ab 



= {tjf\y) U^(y+j)...U; pa (y+(p-l)j) 6(y+pj,x^ 
= (Df)jy\x)\ (5.37) 
Unlike the Laplacian, this is not a Hermitian operator. The source function is now given by 

Rai,(x\y) = F%(x\z) (d ( A (z\w) 5 db 5(w,y). (5.38) 

V / cd 



Since the A matrix may also be written as 



Df) (x\y) = ma Cp(y+(p _ 1) j ) ^^tcac 2(y+ 3 ) m^ (y)5(x5y+ ^ )) 



3 

ab 

T c 2b 



= U a _ c l(y+pj)...Ulf(y+2j) U^(y+j) 5(x,y+pj), 
= U a _J(x)... Ulf(x-(p-l)j)8(x-pj,y), 

D<$) (x\y), (5.39) 

J / ab 
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the source function may be better expressed as 

R ab (x\y) = F$(x\z) (d®) ( z \w) 5 db 5(w,y). (5.40) 

V / cd 

The two key points to note about Eq. (15.401) are that (a) the displacement is applied first 
to the (^-function, and then the smearing is applied, and (b) a displacement in the —k 
direction must be made. 

In summary, to compute 



V(iP, 1>) {D\ p) ^) a Or) (ipD^j (V) ex P Q[U) ip) = G ab (x\y) det Q[U), (5.41) 

one first forms the following source field 

R ab (x\y) = F^{x\z) (£>!?!]) (z\w) 6 db 6(w,y), (5.42) 

V / cd 

= F ac (x\z n ) . . . F C3C2 (z 3 \z 2 ) F C2C1 (z 2 \zx) [D {p) k ) (z^w) 5 db 5(w,y)[5A3) 

; \ / c\d 

solves the linear system of equations 

Qac(x\z) P cb (z\y) = R ab (x\y), (5.44) 
then applies the smearing and displacement at the sink to the result: 

G ab (x\y) = (Df) (x\w)F%(w\z) P cb (z\y), (5.45) 

V J J ad 

= (Df)^ (x\w) F dCni (w\z n( ) . . . F C3C2 (z 3 \z 2 ) F C2Cl (z 2 \ Zl ) P ci6 (^|y),(5.46) 

where the smearing kernel is 

F ab {x\y) = 5 ab 5(x,y) + £ A ab (x\y), (5.47) 

= (l-6^5(x,y)+Z^Uf(x) 5(x+j,y) + U» a (y)* 8(x-j, y)}.(5.48) 

In short, 



at the sink: smear then displace in desired direction, 

at the source: start with point source, displace in opposite direction, then 
smear. 



5.2.2 Chroma walk-through 



Our generalized three-quark propagators were generated using the Chroma QCD library 
developed at JLab [Edw05]. We first used Chroma to make four smeared, 
covariantly-displaced single-quark sources: 
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1. no displacement 

2. +y-direction displacement 

3. — ^-direction displacement 

4. +z-direction displacement 

Each source object was then propagated to all sink sites on the lattice, resulting in a 
smeared, covariantly-displaced source and and unsmeared, undisplaced sink. 

It is in this step that the computationally expensive quark matrix inversion Q~ l is 
performed. This is done using the bi-conjugate gradient algorithm |GL89j . The resulting 
'one-to-all' propagators take the form of Q^q, requiring the inversion of only a single 
column of the quark matrix Q[C/]. 

The sinks were then smeared and covariantly-displaced. For each of the four input 
propagators, there were seven output propagators, one for each of the types of 
displacement: no displacement, ±x, ±y, ±z. These 28 single-quark propagators were then 
combined, three at a time, to make all of the needed generalized three-quark propagators. 
These propagators were saved to disk for later use. 

Combining the indices of the single quark propagators to form the three-quark propagators 
was one of the most computationally intensive steps. This is because the 28 single-quark 
propagators combine into a rich variety of three quark propagators. At the source, there 
can be one of five three-quark elemental operators, one for each type (remembering that we 
have fixed a canonical orientation). At the sink, the number of orientations depends on the 
symmetry of the elemental three-quark operator: 

• single-site (SD): 1 orientation, 

• singly-displaced (SD): 6 orientations, 

• doubly-displaced-I (DDI): 3 orientations, 

• doubly-displaced-L (DDL): 12 orientations, 

• triply-displaced-T (TDT): 12 orientations. 

In addition, we must choose which source quark propagators to each sink quark, and take 
into account that the up and down quarks are indistinguishable (isospin symmetry: 
m u = m d ). 

For the correlation matrix involving a singly-displaced operator at the source (SD) and a 
singly-displaced operator at the sink, there is only one orientation for the quarks at the 
source (by convention), there are six ways to orient the sink quarks, and there are two ways 
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SS 


SD 


DDI 


DDL 


TDT 


ss 


1 


1 


1 


1 


1 


SD 


6 


12 


18 


18 


18 


DDI 


3 


9 


18 


18 


18 


DDL 


12 


36 


72 


72 


72 


TDT 


12 


36 


72 


72 


72 



Table 5.2: The number of three-quark propagators needed for each part of the full nucleon 
correlation matrix. The total number of propagators needed is 671. The diagonal correlation 
matrix elements required only 175 propagators, which were reused for the full run. 



to propagate the displaced source quark: either to the displaced sink quark, or to one of 
the (indistinguishable) undisplaced quarks. 

For example: the number of three-quark propagators required for a full correlation matrix 
run including only single-site and singly-displaced nucleon operators can be calculated as 
follows: 

• SS — SS = 1x1 possibility + 

• SS -~SD = lxl possibility + 

• SD — SS = 1x6 possibilities + 

• SD — SD = 1x6x2 possibilities = 20 three-quark propagators 



We restricted the scope of this project to the iV + (nucleon) channel. Our final correlation 
matrix run used all five types of operators at the source and sink and required 671 
three-quark propagators (see Table I5.2.2j) . These three-quark propagators were stored on 
disk, each file containing one three-quark propagator measured on one configuration. We 
were able to reuse 175 of those three-quark propagators from our diagonal runs. 



5.2.3 Correlation matrix elements 



The C++ program we wrote to tie the three quark propagators together into correlation 
matrix elements was based on the JLab ttt prograuQ. 

The program would first read in the pre-calculated coefficients. For each configuration, or 
bin, the program would read in the three-quark propagator files for that configuration and 
combine them into the desired correlation matrix elements. 



The program name ttt stands for 'tie them together.' 
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We designed the program to process multiple irreps and to generate only the user-selected 
elements of a correlation matrix. Early in the operator tuning and selection process, we 
were only interested in the diagonal elements of the correlation matrices (for effective mass 
functions). Later, we utilized the off-diagonal elements of our correlation matrices to 
extract excited states. 



Chapter 6 



Quark field and gauge link smearing 
parameter tuning 



The smearing method we use to smear our quark fields and gauge links was developed in 
Chapter H] and is summarized here: 



where 



$ = Sip, ip = Sty (6.1) 



a 2 



S(a,n a}Pl n p ) = ^1 + — Aj (6.2) 

&x,V = Yl (Uk( X ) 6 x+'k, V ~ 6 *,v) ( 6 - 3 ) 

fc=±l,±2,±3 

Uk = U [ ^ p] (6.4) 

U [ k n+1 \x) = exp(i P eW(x))uP(x), (6.5) 

Q k (x) = % - (Vi\{x) - n fc (z)) - ^Tr - ft fc (x)) (6.6) 

Qk(x) = C k (x)U k (x) ( no summation over k) (6.7) 

C k (x) = ^(u^Ukix + jpjix + k) (6.8) 

+ U}(x- j)U k (x - j)U 3 (x - j + fc)) . (6.9) 



The four parameters to tune in this smearing procedure are the real (dimensionless) 
Gaussian radius a and staple weight p, and the integer quark field and gauge link smearing 
iterations n a and n p . Due to Chroma conventions, the Gaussian radius values reported 
here will be denoted by a s and related to the a in our equations by 

a s = \plo. 
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6.1 Criterion for judging the effectiveness of smearing 



Recalling our definition of the effective mass function from chapter El a good measure of 
the signal of an operator Oi is the effective mass associated with the correlation function 



given by 

a T M(r) 



= In 



In 



c ii (T) = (n\o i {T)o i (o)\n), 



CJt) 



Cu(t + Or) 





\ c i,l 


2 e - El r + 


Cj,2 


2 e -E 2 r + . . . 


\ c i,l\ 


2 e -E 1 (r+a T ) _|_ 


Cj,2 


2 e -E 2 (T+a T ) -)-... 



i T E x (1 + 0(|q, 2 



2 e -(B 2 _ Bl ) T 



)) 



A reasonable way to proceed is to tune our operators to reduce | c i; fc | 2 / 1 c^^i | 2 for large k, 
because an objective measure of the quality of an operator is the time r at which its 
effective mass function reaches its plateau. In addition, we would like to reduce the 
variance of the Monte Carlo estimate. We define the noise of the effective mass function at 
each time as the absolute ratio of the signal to jackknife error. 

To compare the effectiveness of different values of the quark field smearing parameters, we 
compared the effective mass a T M(r = 4a T ) for each of the three operators at a particular 
temporal separation r = 4a T . Results using 50 quenched configurations on a 12 3 x 48 
anisotropic lattice using the Wilson action with a s ~ 0.1 fm and a s /a T ~ 3.0 are shown in 
Figure 16.31 



6.2 Contamination and noise in extended baryon 
operators 

For single-site (local) hadron operators, it is well known that the use of spatially-smeared 
quark fields is crucial. For extended baryon operators, one expects quark field smearing to 
be equally important, but the relevance and interplay of link-field smearing is less well 
known. Thus, we decided that a systematic study of both quark field and link- variable 
smearing was warranted. 

From our smearing parameter tuning studies [ LMB + 06j we found that Gaussian quark field 
smearing substantially reduces contributions from the short wavelength modes of the 
theory, while stout gauge link smearing significantly reduces the noise from the stochastic 
evaluations. The use of gauge link smearing is shown to be crucial for baryon operators 
constructed of covariantly-displaced quark fields. 
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Figure 6.1: A sample effective mass plot of a single-site operator showing contamination at 
early times due to operator coupling with higher lying modes. 
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Figure 6.2: A sample effective mass plot of a triply-displaced operator showing noise due 
to stochastic gauge link noise. This noise remains present even after significant quark field 
smearing (a s = 4.0, n a = 32). 



Also, the order of approximation of the exponential operator n a determines the maximum 
value of a we can use before the approximation breaks down (see Figure I6T3"]) . Once a 
becomes too large, the smearing operator behaves as (1 + ^-A) n<T — ► ~ cr 2n<T A r, " T , which is 
simply a momentum weighting. Thus, high momentum modes rapidly start to contribute 
even more than in the unsmeared case. 
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6.3 Systematic study of the smearing parameter space 

6.3.1 Trial operators 

Each choice of smearing parameter values necessitates the generation of a new set of 
generalized three-quark propagators. In order to maximize the number of different 
smearing parameter values we could examine, we constructed our three trial operators out 
of three nucleon three-quark propagators having the same operator type at the source and 
sink (single-site, singly-displaced, and triply-displaced-T, see Figure |4~T|) . Our first trial 
operator was a single-site operator Oss projected into the G\ g irreducible representation of 
the cubic point group. Because we did not have different orientations at the sink, we chose 
a particular choice of the Dirac indices for the singly-displaced operator Osd and 
triply-displaced-T operator Otdt- In later runs, we found that the choice of smearing 
parameters was insensitive to the type of extended baryon operator used. 

6.3.2 Quark field smearing 

Without gauge link smearing, the displaced operators were found to be excessively noisy, 
making a meaningful comparison impossible. For this reason, we first tuned the link 
smearing parameters to minimize the noise of the effective mass at the fourth time slice, 
using unsmeared quark fields. We then fixed the gauge link smearing parameters, and 
varied the quark field smearing parameters. These results are shown in Figure 16.31 and 
include gauge-field smearing with n p = 16 and n p p = 2.5. One sees that a T M(r = 4a T ) is 
independent of n a for sufficiently small a s . For each value of n a , a T M(r = 4a r ) first 
decreases as a s is increased, until the approximation to a Gaussian eventually breaks down, 
signaled by a sudden rapid rise in a T M{r = 4a r ). This rapid rise occurs at larger values of 
a s for larger values of n a . 

6.3.3 Gauge link smearing 

Next, we fixed the quark field smearing parameters and studied the effect of changing the 
gauge-field smearing parameters. Before tuning the trial operators, we examined the the 
effective mass a T E(0) associated with the static quark-antiquark potential at a spatial 
separation R = 5a s ~ 0.5 fm and at a particular temporal separation r = to get an idea 
of the effectiveness of different values of p and n p . The results are shown in the leftmost 
plot in Figure 16.41 The behavior is qualitatively similar to that observed when we tuned 
the gauge link smearing parameters for our trial operators, shown in Figure 16 .31 One sees 
that the r = effective mass is independent of the product n p p for sufficiently small values 
of n p p. For each value of n p , a T E(0) decreases as n p p increases, until a minimum is reached 
and a rapid rise occurs. The onset of the rise occurs at larger values of n p p as n p increases. 
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Note that a T E(0) does not decrease appreciably as n p p increases above 2.5. Hence, 
n p p = 2.5 with n p = 16 are our preferred values for the link smearing at lattice spacing 
a s ~ 0.1 fm, based on the static quark-antiquark potential. 

Somewhat surprisingly, we found that changing the link-smearing parameters did not 
appreciably affect the mean values of the effective masses of our three nucleon operators. 
However, the effect on the variances of the effective masses was dramatic. The relative 
jackknife error i](4a T ) of a T M(4a T ) is shown against n p p in the right three plots in 
Figure 16.41 and amazingly, this error shows the same qualitative behavior as in Figure 16.31 
and the leftmost plot in Figure 16.41 One key point learned here is that the preferred 
link-smearing parameters determined from the static quark-antiquark potential produce 
the smallest error in the extended baryon operators. 

The effective masses shown in Figure 16.51 also illustrate these findings. The top row shows 
that applying only quark field smearing to the three selected nucleon operators significantly 
reduces couplings to higher-lying states, but the displaced operators remain excessively 
noisy. The second row illustrates that including only link-field smearing substantially 
reduces the noise, but does not appreciably alter the effective masses themselves. The 
bottom row shows dramatic improvement from reduced couplings to excited states and 
dramatically reduced noise when both quark field and link-field smearing is applied, 
especially for the extended operators. 

6.3.4 Excited state considerations 

Incorporating both quark field and link-variable smearing is crucial for extracting the 
baryon spectrum using gauge-invariant extended three-quark operators. Gaussian quark 
field smearing dramatically diminishes couplings to the short wavelength modes of the 
theory, whereas stout-link smearing drastically reduces the noise in operators with 
displaced quarks. Preferred smearing parameters a s = 4.0, n a = 32, n p p = 2.5, n p = 16 
were found for a lattice spacing a s ~ 0.1 fm and were independent of the baryon operators 
chosen. After examining the behavior of the excited state signal (using the principal 
effective mass method), we decided to reduce the quark smearing in order to decrease the 
noise in the excited states (see Figure [676]) . 

The final values of smearing parameters we selected were: 

• Gaussian smearing: cr s = 3.0, n a = 32 iterations, 

• Stout-link smearing: n p p = 2.5, n p = 16 iterations. 

An issue which remains to be addressed in future work is the dependence of the preferred 
smearing parameters on the quark mass. 
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Figure 6.3: The effective mass a T M(4a r ) for the operators Oss, Osd, Otdt against smear- 
ing radius a s for n a = 1,2,4,8,16,32,64. The gauge field is smeared using n p = 16 and 
n p p = 2.5. Results are based on 50 quenched configurations on a 12 3 x 48 anisotropic lattice 
using the Wilson action with a s ~ 0.1 fm and a s /a T ~ 3.0. The quark mass is such that the 
mass of the pion is approximately 700 MeV. 
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Figure 6.4: Leftmost plot: the effective mass a T E(0) for r = corresponding to the static 
quark- ant iquark potential at spatial separation R = 5a s ~ 0.5 fm against n p p for n p = 
1,2,4,8,16,32. Results are based on 100 configurations on a 16 4 isotropic lattice using 
the Wilson gauge action with f3 = 6.0. Right three plots: the relative jackknife error 
?7(4a T ) of effective masses a T M(4a T ) of the three nucleon operators Oss, @sd, Otdt for 
n a = 32, a s = 4.0 against n p p for n p = 1,2,4,8, 16,32. Results are based on 50 quenched 
configurations on a 12 3 x 48 anisotropic lattice using the Wilson action with a s ~ 0.1 fm, 
a s /a T ~ 3.0. 
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Figure 6.5: Effective masses a T M(r) for unsmeared (black circles) and smeared (red triangles) 
operators Oss, @sd, Otdt- Top row: only quark field smearing n a = 32, a s = 4.0 is used. 
Middle row: only link- variable smearing n p = 16, n p p = 2.5 is applied. Bottom row: 
both quark and link smearing n a = 32, a s = 4.0, n p = 16, n p p = 2.5 are used, dramatically 
improving the signal for all three operators. Results are based on 50 quenched configurations 
on a 12 3 x 48 anisotropic lattice using the Wilson action with a s ~ 0.1 fm, a s /a T ~ 3.0. 
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Figure 6.6: The effects of different values of the quark smearing radius a s on the excited 
states. Throughout, stout-link smearing is used with n p p = 2.5,n p = 16, and n a = 32 
quark smearing interactions are used. The black circles have a s = 4.0, the red squares have 
<j s = 3.0, and the blue triangles have a = 2.0. A Gaussian radius of a s = 3.0 was chosen as 
a balance between high state contamination vs. stochastic noise in the excited states. 



Chapter 7 

Baryon operator pruning 



At this point, we have constructed a set of gauge-invariant single-site and extended baryon 
operators endowed with specific momentum, flavor, and (lattice) spin quantum numbers. 
We have tuned the quark field and gauge link smearing parameters to reduce 
contamination from high-lying modes, and to reduce the noise in the extended operators. 
These efforts have rewarded us with a large set of candidate operators which can be used to 
extract the low-lying baryon spectrum. The numbers of baryon operators of each type 
which project into the different irreducible representations of 0® are given in Table 17.11 

For some set of candidate operators {Oi}, the Hermitian correlation matrix C(r) is the 
fundamental quantity we estimate using the Monte Carlo method: 

Ctf(r) = (01^(^(0)10), 

where we require r > (c.f. Subsection 12.8. ip and once again assume a sufficiently large 
spatial extent T>t and an action which satisfies reflection positivity. The variational 
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Operator type 


Gig 


H 9 


G 2 g 


Gig 


H 9 


G 2 g 


Gig 


Hg 


G 2 g 


Gig 


H 9 


G 2 g 


Single-Site 


1 


2 





4 


3 





3 


1 





4 


1 





Singly-Displaced 


14 


20 


6 


38 


52 


14 


24 


32 


8 


34 


44 


10 


Doubly-Displaced-I 


12 


16 


4 


36 


48 


12 


24 


32 


8 


36 


48 


12 


Doubly-Displaced-L 


32 


64 


32 


96 


192 


96 


64 


128 


64 


96 


192 


96 


Triply-Displaced-T 


32 


64 


32 


96 


192 


96 


64 


128 


64 


96 


192 


96 


Total 


91 


166 


74 


270 


487 


218 


179 


321 


144 


266 


477 


214 



Table 7.1: The numbers of operators of each type which project into each row of the G\ g , 
Hg, and G 2g irreps for the A ++ , S + , N + , and A baryons. The numbers for the H u , and 
G 2u irreps are the same as for the Gi g , H g , and G 2g irreps, respectively. The S° operators 
are obtained from the S + operators by making the flavor exchange u <-> s. The 0~ operators 
are obtained from the A ++ operators by making the flavor replacement u — > s. 
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method for extracting the excited states discussed in Chapter [2] requires us to solve: 

C- 1 / 2 (r )C'(r)C- 1 / 2 (ro)^ = A iUi 

for Ui and Aj. In this study, we will attempt to extract the first eight excited levels in each 
irrep channel. If we were to construct a 179 x 179 correlation matrix for the N + (nucleon) 
G\g channel, then we would have 179 principal correlation functions Aj, of which we would 
only be interested in the first eight. Most of the other 171 principal correlation functions 
would have a negligible signal-to-noise ratio, suggesting that the matrix inversion step 

C(t ) -> C- 1 / 2 ^) 

would most likely fail before we could even solve the eigenvalue problem. Therefore it is 
infeasible to use all of the candidate operators to construct our correlation matrices. We 
chose instead to prune our set of candidate operators down to sixteen operators in each 
irrep before attempting to extract the excited states in those channels using the variational 
method. We feel that this number provides a good balance between the need to use a large 
enough set of operators to reliably extract the states of interest, and the need to use a small 
enough set of operators such that the matrix inversions are not adversely affected by noise. 

Additionally, for each irrep we used the same candidate operators for the even and 
odd-parity channels (Gi g /Gi u , H g /H u , G^g/G^u)- This requirement allows us to improve 
our estimates of the correlation matrix elements by averaging with the backward 
propagating state information in the opposite-parity channel (discussed in 
Subsection I5.1.3p . The two criteria we used to prune our operators were signal quality and 
linear independence. For this work, we focused on the N + (nucleon) channel. 



7.1 Signal quality 

We first calculated the effective mass function associated with each nucleon operator Oj 
using 200 quenched configurations. These effective mass functions required Monte Carlo 
estimates of only the diagonal elements of the nucleon correlation matrices for each irrep: 

• N + Gig/Gi u matrices: 179 diagonal elements, 

• N + H g /H u matrices: 321 diagonal elements, 

• N + Gi g jG2u matrices: 144 diagonal elements. 

The 175 three-quark propagators we generated to construct these elements were saved for 
reuse in our subsequent pruning steps. 

At this pruning stage, we wanted to select operators of each type which had a high quality 
signal. It is reasonable to assume that linear combinations of 'quiet' operators would yield 
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superior results over linear combinations of noisy operators. Examples of quiet vs. noisy 
operators are shown in Figure 17.11 Our measure of signal quality was the effective mass 
function's jackknife error averaged over r = 1, • • • ,16. The consideration of average error 
rather than the signal-to-noise ratio facilitates a direct comparison between the different 
opt 
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Figure 7.1: Effective mass plots for eight representative operators from the complete set of 
extended baryon operators in the G\ u channel. We chose up to ten candidate operators of 
each type (SD, DDI, DDL, TDT) based on the average jackknife error over the first sixteen 
time slices. This process helped us to identify 'quiet' operators (top row, in green) and 
remove some of the noisier operators (bottom row, in red). 



Because there are three single-site type operators in the G\ g jG\ u channels, one in the 
H g /H u channels, and none in the GigjGi u channels, we selected all single-site operators. 
Within each type of extended operator (singly-displaced, doubly-displaced-I, 
doubly-displaced-L, and triply-displaced-T), we sorted the operators in increasing order of 
the average jackknife error over the first sixteen non-zero time separations: r = 1, • ■ ■ ,16 
and selected the first ten operators. At this point, we have selected all of the single-site 
operators and ten extended operators of each type in each irrep. For each operator selected 
in each even-parity channel, we also selected its opposite-parity partner in the odd-parity 
channel, and vice-versa. We performed the signal quality selection each extended operator 
type separately to ensure that we still had a representative set of operators of each type. 
We did not want noise to be the only factor in the selection process; we are primarily 
interested in the ability of an operator to interpolate for the excited states of interest. 
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7.2 Linear independence 

We had eliminated our noisiest operators, but still had over forty candidate operators in 
each channel. We pruned our operators further by looking for subsets of operators which 
excited Hilbert states which were as distinct as possible. We considered such operators 
likely to perform well in the variational method because they would give a more nearly 
orthogonal trial state basis. 

7.2.1 The normalized correlation matrix 

We would like to select a subset of operators such that the corresponding states are as 
distinct as possible. To enable valid comparisons among our candidate operators {Oi}, we 
considered the normalized^ Hermitian correlation matrix elements Cy defined in terms of 

Q. = Cjj(a T ) 

^/Cu(a T )Cjj(a T ) 

(niOie-^Ojin) 
^(niOie-H^OiityiniOje-H^Ojin) 

remembering that the diagonal correlation matrix elements Cu(r) and Cjj(r) are real and 
positive for As in Subsection 12.8.11 we avoid expressions involving the correlation 

matrix evaluated at r = due to the presence of Schwinger terms. 

7.2.2 Pruning by condition number 

For any set of n candidate operators {Oi, O2, • ■ ■ , O n }, we have Ca — 1 V i. In particular: 

Tr((7) = n. 

If all of the states created by our candidate operators were orthogonal, then = Sij, 
which implies that all of the eigenvalues of C are equal to one. In contrast, if two of the 
states were identical, then det((7) = 0, implying that at least one of the eigenvalues of C is 
equal to zero. In order for maintain Tr(C) = n, at least one of the other eigenvalues must 
therefore be greater than one. 

Because our matrix is Hermitian, its eigenvalues are guaranteed to be real. We therefore 
choose as our measure of the degree of linear independence of our operators the condition 



^dn this section we will use a hat to denote the normalized correlation matrix. This is not to be confused 
with the hat notation of Chapter [3] 
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number k, the ratio of the largest to smallest eigenvalue of the normalized correlation 
matrix C: 

^ = / ^max/ / ^min — 1- 

For perfectly orthogonal states k = 1 and for completely degenerate states k — > oo. We 
therefore seek to find a subset of operators in each irrep such that the condition numbers of 
the associated even and odd correlation matrices satisfy some minimization condition. 
Selecting a correlation submatrix having a small condition number has another advantage: 
small condition numbers lead to the rapid convergence of matrix inversion 
algorithms |GL89j . 



7.2.3 Linear independence within each operator type 

We proceeded by calculating all correlation matrix elements having the same type of 
extended operator at the source and sink. We reused the previously generated three-quark 
propagators for the diagonal elements, and generated the additional three-quark 
propagators needed for the off-diagonal elements. 

For each correlation matrix, we calculated the condition numbers of all possible 5x5 
submatrices using the singular value decomposition algorithm |PFTV92] . Because we 
wanted to select the same operators in both the even and odd-parity channels, we 
proceeded as follows for each irrep and operator type separately: 

• Let x denote a subset of five extended operators of a specific type (i.e. all SD or all 
DDL): 

X = { O jj , O j 2 , ■ ■ • , Oi 5 } • 

• Let k 9 (x) and k u {x) be the condition numbers of the even and odd-parity correlation 
submatrices corresponding to the operator subset x, respectively!!. 

• Let n>(x) = max(K 9 (x), k u (x)), the larger of the even and odd-parity condition 
numbers. 

• We selected the operator set x* which minimized k > (x): 

k>(x*) = min(K > (x)) = min (max(K 9 (i), k u (x))) . 



By minimizing the larger of the even and odd-parity condition numbers, we found a 
reasonably distinct set of operators across both the even and odd-parity channels. As a 
check of our results, we examined the smallest eigenvalue A m i n of our selected submatrix to 
verify that it was much greater than its jackknife error, avoiding the possibility of 
signal-to-noise issues. 



2 Remember, we had selected the same operators for both the even and odd-parity channels during our 
signal quality pruning step. 
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A natural question which arises during the application of this method is: how much did we 
lose by requiring that the even and odd-parity channels of each irrep contain the same 
operators? In order to ensure that we were not sacrificing too much, we calculated the 
percent deviation of the condition number from the minimum possible condition number 
for the even and odd-parity channels separately: 

n 9 (x*) — min x ( k 9 (x) ) k u (x* ) — min^, (k u (x)) 
mm x (K 9 (x)) mm x (K u (x)) 

We found that all of these percentages were below 20%, and most were below 5%. The 
extended operator results are shown in Table 17.21 Once again, we selected all single-site 
operators. 



Operator type 


G\g Gi u 


Hg H u 


Gig Glu 


Singly-Displaced 


12.96(21) 12.90(18) 
3% 10% 


1.6972(58) 1.6920(59) 
3% 0% 


5.74(14) 5.46(16) 
0% 0% 


Doubly-Displaced-I 


3.849(30) 3.477(28) 
1% 7% 


1.2956(32) 1.2907(32) 
0% 0% 


2.922(18) 2.842(21) 
1% 0% 


Doubly-Displaced-L 


10.70(16) 10.32(13) 
18% 12% 


4.157(27) 3.894(27) 
16% 0% 


1.697(16) 1.712(17) 
0% 0% 


Triply-Displaced-T 


4.209(42) 4.048(48) 

4% 5% 


1.7796(95) 1.9163(89) 
0% 0% 


1.744(12) 1.856(13) 

4% 2% 



Table 7.2: The condition numbers k and for the 5x5 extended operator submatrices cho- 
sen. The percentages denote how much greater each condition number is than its minimum 
possible value. 



7.2.4 Linear independence across operator types 

In the previous pruning step, we selected all single-site operators and five extended 
operators of each type by considering correlation matrix elements between operators of the 
same type. We next considered correlation matrix elements between operators of different 
types. The total number of three-quark propagators needed was 671. Once again we reused 
previously generated three-quark propagators when possible, and generated the additional 
three-quark propagators needed. 

After normalizing the correlation matrices, we took the best 16 x 16 submatrix according to 
the condition number minimization condition presented in the last subsection. Combining 
the five operators selected from each of the four types with the single-site operators gave 
the following correlation matrix sizes and submatrix numbers for each channel: 

• N + Gig/ Giu matrices: 23 x 23 elements and 245, 157 submatrices, 

• iV + Hg/H u matrices: 21 x 21 elements and 20,349 submatrices, 
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• iV + G 2g /G2u matrices: 20 x 20 elements and 4,845 submatrices, 
where each submatrix has 16 x 16 elements. 

The condition numbers corresponding to the final selected set of sixteen candidate 
operators are presented in Table 17.31 and the new operator number information for our 
final set of sixteen operators is presented in Table 17.41 The operator identification numbers 
used in our programs are given in Appendix 18.61 For completeness, the effective mass plots 
of the operators as calculated from the diagonal correlation matrix elements are shown in 
Figures 17.21 17.31 17.41 17.51 17.61 and 17.71 We also examined the principal effective plots to 
ensure that our selected operators interpolated well for the excited states (see Chapter [8]) . 



Irrep 


K 


(/€ ^min) /^min 


Gig 


68.5(1.5) 


0% 


Giu 


63.85(68) 


4% 


H 9 


53.57(84) 


6% 


H u 


50.09(67) 


0% 


G 2 g 


58.61(76) 


0% 


G 2u 


56.5(2.3) 


3% 



Table 7.3: The condition numbers k for the 16 x 16 final operator submatrices chosen. The 
percentages denote how much greater each condition number is than its minimum possi- 
ble value. The operators in each odd-parity irrep are the opposite-parity partners of the 
operators in the corresponding even-parity irrep. 



Operator type 


Gig/ Glu 


N+ 

Hg/H u 


G 2 g/ G 2u 


Single-Site 


1 








Singly-Displaced 


3 


3 


4 


Doubly-Displaced-I 


5 


3 


3 


Doubly-Displaced-L 


2 


5 


4 


Triply-Displaced-T 


5 


5 


5 


Total 


16 


16 


16 



Table 7.4: The numbers of N + operators of each type which are used in each correlation 
matrix. The operators selected in the G\ u , H u and G 2u irreps are the same as those selected 
for the Gig, H g , and G 2g irreps, respectively. 
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G 1g Effective Masses 
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Figure 7.2: The effective masses for the final sixteen G\ g operators selected by our pruning 
process. 
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H g Effective Masses 
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Figure 7.3: The effective masses for the final sixteen H g operators selected by our pruning 
process. 
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G 2g Effective Masses 
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Figure 7.4: The effective masses for the final sixteen G2 g operators selected by our pruning 
process. 
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G 1u Effective Masses 
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Figure 7.5: The effective masses for the final sixteen G\ u operators selected by our pruning 
process. 
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H u Effective Masses 
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Figure 7.6: The effective masses for the final sixteen H u operators selected by our pruning 
process. 
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G 2u Effective Masses 
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Figure 7.7: The effective masses for the final sixteen Gm operators selected by our pruning 
process. 



Chapter 8 



Baryon operator results: The nucleon 
spectrum 



Now that we have selected sixteen candidate nucleon operators {Oi, O2, ■ ■ ■ , Oiq} in each 
irreducible representation, we conclude this work by extracting the lowest eight energy 
levels in the Gi g , H g , G 2g , Hu, and G 2u nucleon channels, and the lowest seven energy levels 
in the G\ u nucleon channeo 



8.1 Fixed-coefficient correlation functions 

In each irrep channel, we evaluated the elements of a 16 x 16 correlation matrix C(r) 
on 200 quenched configurations. Our next step was to apply the variational method 
discussed in Subsection 12.8.11 to find linear combinations of operators 

16 
o=l 

which couple strongly to single (distinct) states of interest and weakly to the other states. 
In other words, we are looking for fixed coefficient vectors such that: 

16 

(e fe (r)e fc (0)> = <kC ab (r)v bk = vlC(r)v k} (8.1) 

a,b=l 

w \(k\Q k \Q)\ 2 e- EkT , r min <r<r max , 

where [r min , r max ] specifies the range of r values over which the single-exponential 
approximation is valid. It is not guaranteed that such a range exists, but this simple 

1 We were unable to find a fit range which gave a satisfactory estimate for the 8th excited state in the 
Giu channel 
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approach served us well in this exploratory stud)@. A more sophisticated approach would 
attempt to take into account some of the exponential contamination terms neglected by 
this approximation. 

If we solve 

[C-V^CtfC-V^T^UifaTo) = A J (r,r ) Ml (r ) r ), 

for Ui(r, To) and Aj(r, To), where tq is a fixed reference time (we chose To = a T in this work), 
then the principal correlation functions Xi(r, tq) behave as: 



limA fe (r,r ) = e- Ek{T ~ To \l + 0(e 



-A fc (T-T )\ 



A fc = mm \Ei - E k \, 
i^k 

where the E k are the energies of the first n states accessible from the vacuum by the action 
of our trial candidate operators {O a }. The eigenvectors Uiij, r ) may be chosen to be 
orthonormal, in which case we may write the principal correlation functions as: 

A,(r, T ) = u\(r, r ) [C- 1 / 2 (t )C(t)C- 1 / 2 (t q )] u t (r, r ). (8.2) 

Comparing Eqn. 18.21 to Eqn. 18.11 we see that a reasonable choice for the fixed coefficient 
vector Vi is: 

V i = C- 1 '*to)u(T>,To), 

where r* is some fixed time. We found that our estimates for Vi stabilized around r* = 3a T , 
which signaled the reduction of contamination coming from the higher-lying levels. The 
chosen value of r* = 3a T was also small enough to ensure that our estimate for the fixed 
coefficient vectors Vi was based on data having a good signal-to-noise raticH. 

Our fixed- coefficient correlation functions were therefore given by the diagonal elements of 
the rotated correlation matrix, which is denoted by C(t): 

C kk (r) = (0 fc (r)e fc (O)> 
= v\.C{f)v k 

= u\{t = 3a r ,r = a T ) [C^ 1 ^ 2 (a T )C(3a T )C~ 1 ^ 2 (a T )~\ u k (r = 3a T ,r = a T ). 

These fixed-coefficient correlation functions can be used to define fixed-coefficient effective 
mass functions 

C kk (r) 



a T Mk{T) = In 



C kk (r + a T ) 



2 We found that two-exponential fits to our rotated correlation matrix elements v\C(r)vj had inferior 
stability and errors when compared to our single-exponential fit results. The reason for this seems to be our 
effective removal of contamination from the correlation matrix (e.g. see Figure [873]) 

3 The signal-to-noise ratio of Monte Carlo estimates of C(r) decays exponentially with r for baryon 
correlation functions |Lep| 
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8.2 Fitting range 

Our next task was to find the range [r m i n , r max ] (if any) within which the fixed-coefficient 
correlation functions Ckk{ r ) — v \C{ T ) v k behave as single exponentials: 



C kk (r)^\(k\Q k \Q)\ 2 E~ E ^ 



and within which the off-diagonal elements of the rotated correlation matrix 
Cij(r) = v\C{r)vj are consistent with zero: 

C«(r)«0, i^j. 

The vanishing of the off-diagonal elements tells us that any plateaus we observe in our 
effective mass plots (within that range) correspond to orthogonal states. 



8.2.1 Diagonal elements 

In Chapter 0, we discussed the fact that correlation functions will receive contamination 
from the backward propagating states. Because chiral symmetry is spontaneously broken, 
the backward propagating states will have different masses than the forward propagating 
states: 

C£(T) = C^(T-r), but 
C ( k 9 k \r) * cj*>(T-T), 

and similarly for the odd-parity correlation functions. T = a T N T denotes the temporal 
extent of our lattice. 

A sophisticated fitting approach would seek to fit the elements of C{r) and C(T — r) 
simultaneously using independent parameters. The collaboration is currently working to 
develop such a fitting package to use with the production run data. For this work, we 
performed single-exponential fits to the diagonal elements of the rotated correlation matrix, 
neglecting the backward propagating state effects. 

In order to get a rough idea of when the backward propagating state contamination 
becomes significant, we can plot the effective mass function for a trial correlation function 
of the form: 

This simple form assumes the same mass for both the forward and backward propagating 
states, and a coefficient of one for both terms. 
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The associated effective mass is: 

a T M(r) = 



In 
In 



Ctrial(T) 



CtriallT + a 7 



+ e 



-m(T-T) 



g-m(r+a r ) _|_ g— m(T— (r+a T )) 

]_ _|_ e -a T m(N T -2T/a T ) 



a T m + In 



]_ _|_ g— a T m(A r T ~2(T/a T +l) 



5.3) 
U) 
5.5) 



and is plotted in Figure 18 . 1 1 For our lattice, N T = 48 and the lowest level seen in our 
preliminary principal effective mass plots belonged to a G\ g state having a T m « 0.24 
(identified as the proton). From Figure EH] we conclude that in the (worst case) scenario in 
which our backward propagating state is the proton, we will begin to see backward 
propagating state contamination at 



20a T , 



a value which is comparable to the threshold at which our signal-to-noise ratio has 
degraded below acceptable levels for a fit (e.g. see Figure . 
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a T m = 0.24 
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t/a T 

Figure 8.1: The effective mass function for the trial correlation function C(r) = e~ mT + 
e -TO(T-r)_ rj^g numD er of lattice sites in the temporal direction is iV r = 48, and the lowest 
baryon state is taken to be a T m = 0.24. The backward propagating state significantly 
contaminates the effective mass function at times greater than r w 20a r . 



8.2.2 Off-diagonal elements 



The fixed-coefficient correlation functions Cfcfc(r) differ from the principal correlation 
functions Afc(r, r ) because we are only solving the eigenvalue problem once, at r* = 3a T . 
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Thus, we are only guaranteed orthogonality at r* = 3a T : 

The lack of orthogonality for r > r*, no matter how small, means that our fixed-coefficient 
correlation functions will eventually relax to the lowest excited state 

C k k{r) |ci| 2 exp(-E 1 r). 

Consequently, we seek a range of r values over which 

|4-(r)| 2 = Vi/j 

within errors. This would imply that the Vi coefficient vectors are good estimates for the 
coefficients of operators which excite orthogonal states, and would verify that we are not 
excluding too much information by ignoring the off-diagonal elements in our fits. We found 
that this is indeed the case for r values less than r « 30a r . An example of one such 
off-diagonal element is shown in Figure 18.21 

0.01 

0.005 



^ 



-0.005 



"'"'0 10 20 30 40 50 

i/a T 

Figure 8.2: An example of an off-diagonal element of the rotated correlation matrix Cij(r) = 
vjC(r)vj having % = and j = 1. Because we are using fixed coefficients to rotate the matrix, 
we are unable to maintain the orthogonality of the states at all times r. Here we see that we 
can safely work with the fixed-coefficient correlation functions until r « 30a T . After that, 
they all relax to the lowest excited state: C^fr) — > |ci| 2 exp(— E\t). 

Considering the backward propagating state contamination along with the orthogonality of 
the fixed-coefficient correlation functions, we decided to require all of our fit ranges to lie 
within 

"7"min ^rj ^max 20fl r . 

In performing our fits (discussed below), we found that noise was the limiting factor for 
signal quality past r w 17 a r . 
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8.3 Fit method 



To estimate the value of an energy level a T Ek, we fit directly to the fixed-coefficient 
correlation function, not the principal correlation function or any effective mass function^. 
Fitting to correlation functions directly is simpler and introduces less noise into the process. 



8.3.1 Correlated x 2 fitting 

In this study, we performed a single-exponential correlated x 2 fit to the fixed-coefficient 
correlation functions. We fit each fixed-coefficient correlation function Ckk{i~) 
independently. For simplicity in notation, we consider one such fixed-coefficient correlation 
function and denote it by C(r) for the remainder of this subsection. We note that C(t) is a 
function in this subsection, not a matrix. 

The particular form of our fitting function was 

C fit (r, A, a T E &t ) = y4exp(-a T E fit (r/a T )), 

where A and a T E^ are the two fit parameters. This simple form does not take into account 
the periodicity of the lattice; we required our fit ranges to lie within r min = a T and 
Tmax = 20a r to avoid the effects of backward propagating state contamination. 

In our correlated x 2 fit) X 2 is given by 

X 2 (A, a T E m ) = £[C(r) - C flt (r, A, a T E^)]a^ T ,[C{r') - C fit (r', A, a T E at )} (8.6) 

t,t' 

where a' 1 , is the inverse of the covariance matrix: 

1 200 

= N{N _ 1) E^( r ) - C(r)}[C n (r') - C(r% (8.7) 
where C n {r) is the value of C(r) evaluated the n th configuration, and 



1 200 

C(t) = — YCJr). 
K J 200 ^ v ' 

n=l 

The covariance matrix a T)T > measures the amount of correlation between measurements at r 
and r' . 

We used the Levenberg-Marquardt [PFTV92] non-linear fitting algorithm to find the values 
of A and a T Efa which minimized our x 2 ■ To estimate the uncertainty in our fit values, we 



J We denote our energy levels by a T Ek because some may be multi-hadron states. 
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first generated 2048 bootstrap samples of size 200 by sampling randomly with replacement 
from the original set of 200 rotated correlation matrices. We repeated our single 
exponential fit (using the same fit range) on each bootstrap sample, and used the standard 
deviation of the resulting set of fit parameters as our uncertainty We also calculated upper 
and lower limits for our plots using the percentile method [Efr82] . The upper limit of a fit 
parameter is chosen such that 84.1% of the bootstrap sample values lie below it. Similarly, 
the lower limit is chosen corresponding to a percentile of 15.9%. If the fit parameters are 
distributed normally across bootstrap samples, then the upper and lower limits correspond 
to the estimate plus and minus one standard deviation, respectively. We found no 
significant asymmetry in our fit error bars. 

For each fixed-coefficient correlation function separately, we chose our fit range to make the 
X 2 per degree of freedom close to 1.0, and to make the quality factor Q as large as possible 
(> 0.1 was desirable). 

We examined the sensitivity of our fits as we increased or decreased the fit range as a check 
of the quality of our results. We also explored the possibility of adding a second 
exponential term to model some of the high-lying contamination, but those fits performed 
poorly compared to the single exponential method. 

8.3.2 Fitting walk-through 

In summary, the exact method we used to analyze our correlation matrices for each irrep in 

G lg , H g , G 2g , G lu , H u , G 2u is as follows (using 

T o = T min = Or> r * = 3a T , T max = 20a r , T = a T N T , and N T = 48): 

1. Average the correlation matrix elements over all 200 configurations to form the best 
estimates for the elements of C(r) using the improved estimation method discussed 
in Chapter [3j 

200 

<%\t) = ^ E i [ C %n(r) + C^:\T - r) + Cg;(r) + C^>(T - r)_ , 

n=l 

where Cj£ (t) is a correlation matrix element measured on the n th configuration 
generated by the Markov chain (see Chapter [3]), and p = 0, 1 for even and 
odd-parity, respectively. 

2. Solve C~ 1/2 (tq)C(t*)C~ iI2 (t q ) u t = Aj u» for the eigenvector ti;. 

3. Calculate Vi = C _1 / 2 (r )nj. 

4. For each configuration n — 1, • • • , 200, store the elements of the rotated correlation 
matrix C for all r values: 

16 

Cij, n {f) = ^ V* a C abtn (T) Vbj . 
a,b=l 
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5. Calculate the mean and jackknife error of the off-diagonal elements of the rotated 
correlation matrix C{r) and verify that they are consistent with zero for 

6. Perform single-exponential fits to each fixed-coefficient correlation function C^r). 
In these fits, the covariance matrix is calculated using Eqn. 18.71 on the sample of 
rotated correlation matrices and then inverted for use in the expression for x 2 given 
in Eqn. 18.61 The fit parameters are varied until \ 2 is minimized. In each fit, choose a 
fit range such that the \ 2 P er degree of freedom is near 1.0, the quality factor Q is as 
large as possible (> 0.1 is desirable), and the fit parameter a T E k ^ t is stable under 
small variations in the fit range. Store the value of a T E k fit as the estimate of the level 
value. 

7. To estimate the uncertainty in each level estimate, generate a large number (we used 
2048) of bootstrap samples of size 200 by sampling randomly with replacement from 
the original set of 200 rotated correlation matrices. Perform a single-exponential fit 
using the same fit range used for the fit on the original sample. Once again, the 
covariance matrix is estimated using Eqn. 18. 7[ this time on the members of the 
bootstrap sample. Each bootstrap sample yields an estimate of the level value a T Ek^ t . 

8. Store the standard deviation of the bootstrap estimates of a T E k ^ t as the uncertainty 
in the level estimate. 

9. Examine the fixed-coefficient effective masses, the level estimates, and the principal 
effective masses to verify that the fits are reasonable and not adversely affected by 
outlying data. 



8.4 Lattice nucleon spectrum results 

This lattice study was performed using 200 quenched configurations on a 12 3 x 48 
anisotropic lattice using the Wilson action with a s ~ 0.1 fm and a s /a T ~ 3.0. The quark 
mass was such that the mass of the pion was approximately 700 MeV. We were able to 
perform satisfactory fits yielding the lowest-lying eight states in each lattice spin-parity 
channel, with the exception of the highest (eighth) G\ u level. Our spectrum fit results are 
tabulated in Tables EH and EM and plotted in Figures EM EM EM EH EH and EH The 
nucleon spectrum levels extracted from these fits are illustrated in Figure 18.91 
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Table 8.1: The final spectrum results for the even-parity channels. These results are based 
on 200 quenched configurations on a 12 3 x 48 anisotropic lattice using the Wilson action 
with a s ~ 0.1 fm and a s /a T ~ 3.0. Our pion mass for this study was a T M n = 0.1125(26) 
(see Figure IzTTl) . or approximately 700 MeV. 
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Table 8.2: The final spectrum results for the odd-parity channels. These results are based 
on 200 quenched configurations on a 12 3 x 48 anisotropic lattice using the Wilson action 
with a s ~ 0.1 fm and a s /a T ~ 3.0. Our pion mass for this study was a T M n = 0.1125(26) 
(see Figure [27TT) . or approximately 700 MeV. We did not find a satisfactory fit range for the 
eighth level in the G± u channel. 
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G 1g Nucleon Fit Results 
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Figure 8.3: The lowest eight levels for the G\ g nucleon channel. The green circles are the 
fixed coefficient effective masses, and the red squares are the principal effective masses. The 
fits were made to the fixed coefficient correlation functions and are denoted by the blue lines. 
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H g Nucleon Fit Results 
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Figure 8.4: The lowest eight levels for the H g nucleon channel. The green circles are the 
fixed coefficient effective masses, and the red squares are the principal effective masses. The 
fits were made to the fixed coefficient correlation functions and are denoted by the blue lines. 
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G 2g Nucleon Fit Results 



03 



1.2 



0.9 



0.6 



0.3 



Level 1 

a T E fi =0.393(16) 




Q ^ /(d.o.f.)=0.71., Q=0-65 



Level 2 

a T E fi =0.409(13) 



7(d.o.f.)=0.63, Q=0.74 



CO 



CO 



CO 



1.2 

0.9 

0.6 

0.3 


1.2 

0.9 

0.6 

0.3 


1.2 

0.9 

0.6 



Level 3 

a T E fi =0.420(13) 



a 



7(d.o.f.)=1.04 f Q=04(f 



Level 4 

a T E fi =0.425(11) 



7(d.o.f.)=0.99, Q=0.43 



Level 5 

. a x E fi =0.586(^3) 



-II 



X 2 /(d.o.f.)=1.10. Q=f.3(i 

L_i i i i I i i i i 1 i i_ 



Level 6 



• a x E fi =0.602(^) 



m 



□ 



o 



% 7(d.o.f.)=0.66. Q=ft-7f 

L_l I I I I 1 I I I ] I L 



0.3 



Level 7 

a T E fi =0.613(|6) 



IP" 
[] □ 



— 



0<) [] 



^c 2 /(d.o.f.)=0.47 f Q=0.89 



Level £ 

a T E fi =0.63i((l$) 



jC 7(d.o.f.)=0.95 ? gHjffi 



5 10 15 

x / a 



5 10 15 

x / a_ 



Figure 8.5: The lowest eight levels for the G23 nucleon channel. The green circles are the 
fixed coefficient effective masses, and the red squares are the principal effective masses. The 
fits were made to the fixed coefficient correlation functions and are denoted by the blue lines. 
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G 1u Nucleon Fit Results 
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Figure 8.6: The lowest eight levels for the G\ u nucleon channel. The green circles are the 
fixed coefficient effective masses, and the red squares are the principal effective masses. The 
fits were made to the fixed coefficient correlation functions and are denoted by the blue lines. 



143 



H u Nucleon Fit Results 
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Figure 8.7: The lowest eight levels for the H u nucleon channel. The green circles are the 
fixed coefficient effective masses, and the red squares are the principal effective masses. The 
fits were made to the fixed coefficient correlation functions and are denoted by the blue lines. 
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G 2u Nucleon Fit Results 
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Figure 8.8: The lowest eight levels for the G2 U nucleon channel. The green circles are the 
fixed coefficient effective masses, and the red squares are the principal effective masses. The 
fits were made to the fixed coefficient correlation functions and are denoted by the blue lines. 
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Nucleon Mass Spectrum 
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Figure 8.9: The low-lying 1=1/2, I 3 =+l/2 nucleon spectrum extracted from 200 quenched 
configurations on a 12 3 x 48 anisotropic lattice using the Wilson action with a s ~ 0.1 fm 
and a s /a T ~ 3.0. The vertical height of each box indicates the statistical uncertainty in that 
estimate. Our pion mass for this study was = 0.1125(26), or approximately 700 MeV. 

A different color was chosen for each level in a symmetry channel to help the reader discern 
among the different levels. 
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Nucleon Mass Spectrum (Exp) 



2.8- 



!_K 113 (2700) K 113 (2700) K 1 13 (2700) 

(x2) (x2) I ^(2600) L ^(2600) L ^(2600) 



2.4 - 



G 17 (2190) G 17 (2190) G 17 (2190) 
P ii( 210 ° ) F is( 2 000) F 15( 200 °) Sl1 (2090) D "^o7o) 
F 17 (1990) F, 7 ( 1990) F 17 (1990) 
P 13 (1900) 



5 1-6 



1.2 - 



0.8- 



0.4- 







H 19 (2220) H 19 (2220) (x2) 



(x2) 

G.J2250) G 19 (2250) (x2) 



15 



(2200) D (2200) 



P 1 ,(1720) 



S.J1650) 



D 13 (1700) 



£ii (1710) F 15 (1680) F 15 (1680) D 15 (1675) D 15 (1675) 

s"7l535) D^fl520) 



P 11 (1440) 



P 11 (939) 



G 



ig 



H 



g 



G 



2g 



G 



1u 



H 



u 



G 



2u 



Figure 8.10: The 7=1/2, J 3 =+l/2 nucleon spectrum as determined by experiment [Y + 06] 
and projected into the space of lattice spin-parity states. Black denotes a four-star state, 
blue denotes a three-star state, tan denotes a two-star state, and gray denotes a one-star 
state. 
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Nucleon Mass Spectrum (RQM) 
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Figure 8.11: The low-lying 7=1/2, I 3 =+l/2 nucleon spectrum up to 2820 MeV as predicted 
by the relativistic quark model |CR93j and projected into the space of lattice spin-parity 
states. 'Missing resonances' are displayed in red, and observed resonances are labeled by 
their assigned state in the experimental spectrum. Note that the location of a labeled state 
is given by its predicted energy value, and the text of the label tells its experimentally 
measured energy value. For a labeled state black denotes a four-star state, blue denotes a 
three-star state, tan denotes a two-star state, and gray denotes a one-star state. 
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8.5 Discussion 

Our final nucleon spectrum results for this study are tabulated in Tables 18.11 and 18.21 and 
illustrated in Figure 18.91 Comparison of these results with experiment is not justified 
because the quenched approximation was used, an unphysically large u, d quark mass was 
used, and the lattice volume was too small. Nevertheless, such a comparison is still 
interesting, so we discuss our results with an eye towards the experimentally observed 
spectrum, shown in Table 18.31 and Figure 18.101 It is also interesting to compare our results 
to the predictions of the relativistic constituent quark model |CI86t ICR93j ICR94j . shown in 
Tables El and Eland Figure EIU 

We emphasize that most of the levels shown in Figure 18.91 have never before been 
calculated using first-principles QCD. Figure 18.91 is actually a first glimpse of the low-lying 
nucleon spectrum as predicted by QCD (with the caveats mentioned above). 

First, the lightest nucleon state has a mass value around 0.24 a~ l and resides by itself in the 
G\g channel, indicating a spin-| state. This clearly corresponds to the proton at 939 MeV. 

Second, there is a cluster or band of odd-parity states just below 0.35 a" 1 and well 
separated from higher-lying states. In this band, there are two G\ u states, three H u states, 
and one Go u state. The G 2u level most likely corresponds to the _Di 5 (1675) spin-| 
resonancqj. This means that one of the H u levels must correspond to this same state (see 
Table I&31) . The other two H u states most likely correspond to the spin-| _Di3(1520) and 
_Di 3 (1700) resonances. The two G lu levels seem to correspond to the spin-^ S , n(1535) and 
Sii(1650) resonances. There are no other experimentally-observed odd-parity nucleon 
resonances below 2.0 GeV. So it appears that there is a one-to-one matching of the states 
in this band with experiment: every state that should be seen is seen, and there are no 
extra states. However, given that the mass of the pseudoscalar pion is about 0.11 a^ 1 , we 
expect a two-particle N-ir 5-wave state in the odd-parity G\ u channel near this energy. We 
are either missing such a state, or are missing one of the S , n(1535) or 5 , n(1650) states. The 
most likely explanation is that the two-particle state is absent from our spectrum due to 
the quenched approximation combined with our use of single particle operators. Future 
unquenched calculations should resolve this question. 

Third, we observe another band or cluster of states in the even-parity channels around 
0.4 a" 1 . There are four states in each of the G\ g , H g , G2 g channels, although two additional 
levels in H g are only slightly higher. The most striking feature of this band is the lack of 
one Gig state lying just below the previously discussed odd-parity band. Experiments 
clearly show a single resonance near 1440 MeV, known as the Roper resonance. The Roper 
resonance in our spectrum occurs at much too high of an energy compared with 
experiment. The most likely explanation for this is our use of unphysically large quark 

Experimentalists refer to a baryon resonance having isospin / and total angular momentum J using 
the notation L21 2J, where L = S, P, D, F,G, H, I, K ■ ■ ■ is the orbital angular momentum of an Nir system 
having the same J p as the state. 
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masses. There is some evidence |MCD + 05] from previous quenched lattice calculations that 
the Roper mass dramatically decreases when the u, d quark mass is decreased such that the 
pion mass falls below 300 MeV (our pion mass is around 700 MeV). 

In the observed spectrum (Figure [8 . 101) between 1.6 GeV and 1.8 GeV, there are currently 
three well-established even-parity resonances: the Fi 5 (1680), the Pn(1710), and the 
Pi3(1720). A second band of four more even-parity resonances has been tentatively 
identified near 2.0 GeV containing the P 13 (1900), the P 17 (1990), the Pi 5 (2000), and the 
Pn(2100) states. Our results do not reproduce this two-band pattern. Rather, we find a 
single band of four states in the each of the Gi g , Hq, and G^g channels around 0.4 a~ x . 

One possibility is that our calculations have identified resonances in the 1.6 to 1.8 GeV 
range which have not yet been observed in experiments. These states would correspond to 
some of the 'missing resonances' predicted by the relativistic constituent quark model (see 
below and Figure IHTTT]) . 

However, if the mass splittings among the experimentally observed levels are particularly 
sensitive to the u, d quark mass, one may speculate that all of those levels, including the 
Roper, might merge to a single band at large u, d mass. If so, one would observe four 
nearly degenerate levels in the G\ g channel, five levels in the H g channel, and three in the 
G2 g channel. Such a pattern would more closely resemble the results of our calculation (see 
Figure EM- 

Additionally, the single-band structure observed in the even-parity channels may split 
when we add a clover (a ■ F) improvement terrrj^l to our action [SW85[ IEHK98b| IEHR03j . 
Such a term not only compensates for the 0(a) discretization errors introduced by the 
chiral-symmetry-breaking Wilson term, but may also restore level splitting behavior 
disrupted by the Wilson term at 0(a). 

Increased statistics and a finite lattice spacing study (involving data from several different 
lattice spacings a s ) are needed in order to deduce the spin content in this even-parity band. 
Note that we also expect two-particle 7r-S' 11 (1535) and 7r-P 13 (1520) states to occur in the 
even-parity spectrum just slightly above this band's energy, but the quenched 
approximation might again make such states inaccessible to the operators used in this 
study. Whether or not all of these levels merge into a single band at large u, d mass will be 
eventually answered; in future calculations we will determine the mass spectrum for various 
values of the u, d quark mass to resolve this issue. 

Although our odd-parity results agree qualitatively with the quark model [CR93J, the 
quark model predicts more even-parity states around 2 GeV than we found. Quark model 
calculations [OR93] predict a rather dense spectrum of four | + states, five | + states, three 
| + states, and one | + state in the 1.65 to 2.15 GeV energy region. If this were correct, 
then our lattice spectrum would contain a low-lying band of five states in the G\ g irrep, 
nine states in the H g irrep, and four states in the G2 g irrep (see Figure IHTTj) . The number 

6 The a matrices are given by — —§[7^,7^], and F is a discretization of the gauge field strength. 
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State and 




Experimental 




Appearance 


in OP 

h 


Irrep 




Mass (MeV) 


J p 


Status 






G2g 


G\u 


H u 


G*2n 


£i(939) 


1 + 
2 




1 

















£11 (1440) 


1 + 
2 


* * 


1 

















£i 3 (1520) 


3- 
2 


* * ** 














1 





Su(1535) 


1- 
2 


* * ** 











1 








5u(1650) 


1- 
2 


* * ** 











1 








£15 (1675) 


5- 
2 


* * 














1 


1 


F 15 (1680) 


5 + 
2 







1 


1 











^13(1700) 


3- 
2 


* * * 














1 





£n(1710) 


1 + 
2 


* * * 


1 

















Pi 3 (1720) 


3 + 
2 


* * ** 





1 














Pi 3 (1900) 


3 + 
2 


** 





1 














F 17 (1990) 


7 + 
2 


** 


1 


1 


1 











Fi5(2000) 


5 + 
2 


** 





1 


1 











£i 3 (2080) 


3- 
2 


** 














1 





5n(2090) 


1- 
2 


* 











1 








£n(2100) 


1 + 
2 


* 


1 

















G 17 (2190) 


7- 
2 


* * ** 











1 


1 


1 


L> 15 (2200) 


5- 
2 


** 














1 


1 


#19 (2220) 


9 + 
2 


* * ** 


1 


2 














Gi 9 (2250) 


9- 
2 













1 


2 





£,n(2600) 


11- 
2 


* * * 











1 


2 


1 


£i,i3(2700) 


13 + 
2 


** 


1 


2 


2 












Table 8.3: The current experimental values [Y + 06] for all known nucleon resonances. The 
state names are given in spectroscopic notation: L 2 j 2j, where L = S, P, D, F, ■ ■ ■ is the 
orbital angular momentum of an Ntc system having the same J p as the state, I is the 
isospin, and J is the total angular momentum. A four-star experimental status implies 
that existence is certain and that the properties are fairly well explored. A three-star status 
implies that existence ranges from very likely to certain, but further confirmation is desirable. 
Two stars implies that the evidence for existence is only fair, and one star implies that the 
evidence is poor. The number of times each state is expected to appear in each lattice OP 
irrep (obtained from subduction) is also shown. The experimental uncertainties are at the 
5% level or less. 
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Model 




Ntt State 


Experimental 




Appearance 


in Of 


Irrep 


State 


J p 


Assignment 


Status 


Gi 9 


H 9 


G 2 g 


G\u 


H u 




5^(1460) 


i- 

2 


5n(1535) 


* * ** 











1 








D\ 3 (U95) 


3- 
2 


£>i 3 (1520) 


* * ** 














1 





5^(1535) 


1- 
2 


5n(1650) 


* * ** 











1 








Pn(1540) 


1 + 
2 


Pn(1440) 




1 

















L>? 3 (1625) 


3~ 
2 


^13(1700) 


* * * 














1 





1^(1630) 


5~ 
2 


£>i 5 (1675) 
















1 


1 


F£(1770) 


5 + 
2 


Pi 5 (1680) 







1 


1 











^(1770) 


1 + 

2 


Pn(1710) 


* * * 


1 

















Pi(1795) 


3 + 
2 


P 13 (1720) 







1 














A 2 3 (1870) 


3 + 
2 


— 


— 





1 














Pi 4 ! (1880) 


1 + 
2 


— 


— 


1 

















Pi 3 3 (1910) 


3 + 
2 


— 


— 





1 














5^(1945) 


1- 
2 


511(2090) 


* 











1 








Pi 4 3 (1950) 


3 + 
2 


— 


— 





1 














.D? 3 (1960) 


3~ 
2 


Pi 3 (2080) 


** 














1 





Pfi (1975) 


1 + 
2 


— 


— 


1 

















^(1980) 


5 + 
2 


F 15 (2000) 


** 





1 


1 











Ff 5 (1995) 


5 + 
2 


— 


— 





1 


1 











P 1 1 7 (2000) 


7 + 
2 


Pi 7 (1990) 


** 


1 


1 


1 











P 1 5 3 (2030) 


3 + 
2 


— 


— 





1 














5^(2030) 


1- 
2 


— 


— 











1 








D 4 3 (2055) 


3 _ 
2 


— 


— 














1 





P 1 6 1 (2065) 


1 + 

2 


Pn(2100) 


* 


1 

















5^(2070) 


1- 
2 


— 


— 











1 








r~»2 /onon\ 

iJ{ 5 (zU8(J) 


5~ 
2 


















1 


1 


G] 7 (2090) 


7~ 
2 


Gi 7 (2190) 













1 


1 


1 


Pf 3 (2095) 


3~ 
2 


















1 





P? 5 (2095) 


5~ 
2 


Pi 5 (2200) 


** 














1 


1 



Table 8.4: The quark model predictions [CR93] for the excited nucleon spectrum below 
2100 MeV. The proton is used to set the parameters of the model. The state names are 
given in spectroscopic notation: L 2 i 2j, where L = S, P, D, F, ■ ■ ■ is the orbital angular 
momentum of an Nir system having the same J p as the state, I is the isospin, and J is 
the total angular momentum. The superscripted integer on the model state name denotes 
the principal quantum number in the quark model (see [CR93] ). The number of times each 
state is expected to appear in each lattice irrep is also shown. Dashes indicate 'missing 
resonances.' 
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Model 




A^7r State 


Experimental 




Appearance 


in 0° 


Irrep 




State 


T P 


Assignment 


Status 


Gig 


H 9 




G\ u 


H u 




(2145) 


I — 

2 















1 








£>«,(2165) 




2 


















1 





D[ 3 (2180) 




2 


















1 





1^(2180) 


c — 


2 


















1 


1 


S' 1 7 1 (2195) 


1 
2 















1 








G? 7 (2205) 


7 

2 















1 


1 


1 


P 1 7 1 (2210) 


2 






1 

















GJ 9 (2215) 


2 


Gi 9 (2250) 













1 


2 





^15(2235) 


c; — 


2 


















1 


1 


G? 7 (2255) 


1 — 
/ 

2 















1 


1 


1 


D? 5 (2260) 


c: — 
O 

2 


















1 


1 


^5(2295) 




2 


















1 


1 


D? 5 (2305) 


c — 


2 


— 


— 














1 


1 


Gf 7 (2305) 


2 















1 


1 


1 


^9(2345) 


Q + 
2 


#19 (2220) 




1 


2 














Gf 7 (2355) 


I 
2 















1 


1 


1 


^(2390) 


7 + 

2 






1 


1 


1 











Ff 7 (2410) 


7 + 
/ 1 

2 






1 


1 


1 











^(2455) 


7 + 
2 






1 


1 


1 











^1,11 (2490) 


11 + 

T 






1 


2 


1 











#£,(2490) 


q + 

y 1 

2 






1 


2 














#i 9 (2500) 


9+ 
2 






1 


2 














H f n (2600) 


11 + 

T 






1 


2 


1 











J 1 1 ill (2600) 


xi — 
T 


Ji,n(2600) 


Sf.Sf.-lf. 





n 
u 


n 
u 


1 


2 


1 


J? 1.(2670) 


H- 
2 















1 


2 


1 


A 3 11 (2700) 


U - 
2 















1 


2 


1 


^,13(2715) 


13- 
2 















1 


2 


2 


A 4 11 (2770) 


11- 
2 















1 


2 


1 


^,13(2820) 


13 + 
2 


^1,13(2700) 


** 


1 


2 


2 












Table 8.5: The quark model predictions [CR93] for the excited nucleon spectrum from 2100 
MeV to 2820 MeV. Dashes indicate 'missing resonances.' 
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of states appearing in the lowest H g band is too small to support this. Future studies will 
extract more levels in the H g channel to shed more light on this issue. 

With statistics based on only 200 configurations, the interpretation of the states above 
0.45 a" 1 is somewhat problematic. However, the absence of states between 0.45 a" 1 and 
0.55 a^ 1 in the G*2g channel is an interesting feature of the spectrum and does appear to 
tentatively agree with experiment. The quark model does not predict such a large gap (see 
Figure 15. lip . On the other hand, our calculation also shows a large clustering of levels 
around 0.5 a~ l in the channel, which agrees more with the predictions of the quark 
model than with the experimental data. 
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8.6 Conclusion and outlook 

In this work, we have presented a systematic approach to baryon operator design for use in 
lattice QCD baryon spectrum calculations. We described the formalism for extracting 
excited states from correlation matrix elements between quantum operators, and 
emphasized the importance of using operators which couple strongly to the states of 
interest and weakly to the high-lying contaminating modes. 

We discussed the evaluation of correlation matrix elements using the Monte Carlo method 
on an anisotropic space-time lattice. The anisotropy allows both a fine temporal resolution 
for the examination of exponentially decaying correlation functions, and a coarse spatial 
resolution enabling the use of a larger spatial volume. The signal quality associated with 
the diagonal correlation matrix elements was defined in terms of the effective mass 
function, which provided not only a quantitative measure of the contamination and noise, 
but also a method of visualizing the presence of stationary energy states as plateaus in 
effective mass plots. 

The baryon operators used in this study were composed of gauge-covariantly-displaced 
three-quark operators. The covariant quark displacements allowed us to incorporate radial 
structure resulting in operators which better interpolated for excited baryon states. We 
found it essential to use smeared gauge links in our extended operators to suppress the 
noise introduced by using gauge links to covariantly displace the quarks. Additionally, we 
found that while quark field smearing had little effect on the noise, it dramatically reduced 
the coupling to high-lying contaminating modes. By systematically exploring the smearing 
parameter space while monitoring the behavior of the effective mass plots associated with 
three trial operators, we found values of the smearing parameters which reduced the noise 
and contamination in the correlation matrix elements while preserving the integrity of the 
excited state signals. 

We then used group theory to construct operators which transformed irreducibly under the 
symmetry group of the cubic spatial lattice. This resulted in operators with well-defined 
lattice quantum numbers, allowing us to identify the J p quantum numbers of the 
corresponding continuum states. The good quantum numbers on our lattice were the linear 
momentum (chosen to be zero for our spectrum calculations), color (our operators were 
gauge-invariant, or colorless), isospin (e.g. nucleon operators only excited nucleon states), 
and the lattice spin-parity irreducible representation label corresponding to the spinorial 
representations of Of, the double- valued octahedral point group of rotations and reflections 
on the cubic spatial lattice. This last label, which took the values G\ g , H g , G\ u , H u , and 
G 2u , provided a way to identify the continuum J p quantum numbers of the spectral states. 

The group-theoretical approach led to an unmanageably large number of extended baryon 
operators. We developed a method to select sixteen operators in each lattice spin-parity 
channel to be used for the final extraction of the spectrum. This pruning method was 
based on three crucial requirements: low noise, maximal independence, and good overlap 
with the excited states as judged by principal effective mass quality. We evaluated the 
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diagonal elements of the operator correlation matrix to remove noisy operators, and then 
selected the sixteen operators whose renormalized correlation matrix at a fixed small time 
separation had a low condition number for both the even- and odd-parity channels. 

We then applied a variational method to rotate our basis set of operators into a new set 
which could be used to extract the lowest seven or eight levels of the spectrum in each 
channel. We found a range of r values over which the temporal correlation function for 
each of these new operators could be fit to a single exponential. After performing the fits, 
we interpreted the results in terms of both experimental data and quark model predictions. 

Although comparison with experiment is not justified, the pattern of levels obtained 
qualitatively agrees with the observed spectrum. We also compared our spectrum to 
relativistic quark model predictions; the quark model predicts more low-lying even-parity 
states than this study found, but both the quark model and this study predict more 
odd-parity states near 2 GeV than currently observed in experiments. These results are not 
only interesting, but also serve as a powerful validation of the operator design methodology 
described in this work and provide a preview of what we expect to achieve in this 
long-term effort. 

We focused on the nucleon channel in this work, but the methodology used here is readily 
applied to the other baryon channels (A, S, 5, A, and Q). Our approach to building 
operators and extracting the spectrum is also expected to be equally effective when 
adapted for the construction of meson operators, an important next step. 

An important improvement to our calculation technique will be the use of all-to-all 
propagators |FJO + 05| . rather than one-to-all propagators, when evaluating correlation 
matrix elements. All-to-all propagators will allow us to treat the source and sink operators 
separately, an improvement over the current method in which we are forced to express 
correlation matrix elements in terms of thousands of three-quark propagator terms. It will 
be computationally feasible to spatially average over the source as well as the sink to 
improve our statistics, and to evaluate correlation matrix elements between multi-hadron 
operators. 

This research was carried out as part of the ongoing Lattice Hadron Physics Collaboration 
QCD spectroscopy project. The LHPC has recently been awarded over ten million CPU 
hours on the QCDOC at Brookhaven National Laboratory, a massively parallel ASIC 
supercomputer [BCC + 04] . We are currently generating 1,000 large- volume unquenched 
gauge configurations. Because this unquenched run includes the fermion determinant, we 
are correctly including quark loop effects in the calculation of our propagators. Using 
clover- improved Wilson fermions, an improved gauge action, and a pion mass around 
300 MeV, we expect this production run data to resolve many of the questions raised by 
the preliminary spectrum results presented here, and to provide us with a first ab initio 
look into the baryon spectrum as predicted by QCD. 



7 QCDOC stands for "QCD On a Chip" and ASIC stands for "Application-Specific Integrated Circuit. 



Appendix: Final operator selection 



The number of nucleoli operators of each type is shown in Table [61 For future studies, we 
record here the corresponding identification numbers for the final sixteen nucleon operators 
selected from each irreducible representation in Tables [3, [HI and [HI These identification 
numbers are used to index the different operators in our projection_coef f icients data 
files. 



iV + Operator type 


Gig 


H g 


G 2 g 


Single-Site 


3 


1 





Singly-Displaced 


24 


32 


8 


Doubly-Displaced-I 


24 


32 


8 


Doubly-Displaced-L 


64 


128 


64 


Triply-Displaced-T 


64 


128 


64 


Total 


179 


321 


144 



Table 6: The numbers of operators of each type which project into each row of the Gi g , H g , 
and Gig irreps for the N + baryons. The numbers for the Gi u ,H u , and G 2u irreps are the 
same as for the Gi g ,H g , and G 2g irreps, respectively. 
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Gig/Gi u Operators 



Operator Number 


Operator Type 


Operator ID 


1 




Single-Site 


2 


2 




Singly-Displaced 


11 


3 




Singly-Displaced 


17 


4 




Singly-Displaced 


20 


5 




T^v 1 1 T~\ 1 IT 

Doubly-Displaced-1 





6 




Doubly-Displaced-I 


4 


7 




Doubly-Displaced-I 


5 


8 




Doubly-Displaced-I 


9 


9 




Doubly-Displaced-I 


12 


10 




Doubly-Displaced-L 


4 


11 




Doubly-Displaced-L 


10 


12 




Triply-Displaced-T 


3 


13 




Triply-Displaced-T 


5 


14 




Triply-Displaced-T 


9 


15 




Triply-Displaced-T 


11 


16 




Triply-Displaced-T 


25 



Table 7: The identification numbers for the final sixteen nucleon operators selected from the 
G\ g jG\ u channels. The ID number corresponds to the operator number within each type 
(see Table E]) as indexed in our projection_coef f icients data files. 
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H g /H u Operators 



Operator Number 


Operator Type 


Operator ID 


1 


Singly-Displaced 


9 


2 


Singly-Displaced 


10 


3 


Singly-Displaced 


31 


4 


Doubly-Displaced-I 


8 


5 


T^v 1 1 T~\ 1 IT 

Doubly-Displaced-1 


1 1 


6 


Doubly-Displaced-I 


31 


7 


Doubly-Displaced-L 


47 


8 


Doubly-Displaced-L 


54 


9 


Doubly-Displaced-L 


84 


10 


Doubly-Displaced-L 


113 


11 


Doubly-Displaced-L 


124 


12 


Triply-Displaced-T 


35 


13 


Triply-Displaced-T 


71 


14 


Triply-Displaced-T 


86 


15 


Triply-Displaced-T 


95 


16 


Triply-Displaced-T 


104 



Table 8: The identification numbers for the final sixteen nucleon operators selected from the 
H g /H u channels. The ID number corresponds to the operator number within each type (see 
Table E]) as indexed in our projection_coef f icients data files. 
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GiglG^u Operators 


Operator Number 


Operator Type 


Operator ID 


1 


Singly-Displaced 





2 


Singly-Displaced 


1 


3 


Singly-Displaced 


2 


4 


Singly-Displaced 


6 


5 


Doubly-Displaced-I 


5 


6 


Doubly-Displaced-I 


6 


7 


Doubly-Displaced-I 


7 


8 


Doubly-Displaced-L 


32 


9 


Doubly-Displaced-L 


37 


10 


Doubly-Displaced-L 


41 


11 


Doubly-Displaced-L 


52 


12 


Triply-Displaced-T 


1 


13 


Triply-Displaced-T 


33 


14 


Triply-Displaced-T 


45 


15 


Triply-Displaced-T 


51 


16 


Triply-Displaced-T 


61 



Table 9: The identification numbers for the final sixteen nucleon operators selected from the 
G2 9 /G2u channels. The ID number corresponds to the operator number within each type 
(see Table E]) as indexed in our projection_coef f icients data files. 
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